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Abstract 

Particles have tremendous potential as astronomical messengers, and conversely, studying 
the universe as a whole also teaches us about particle physics. This thesis encompasses 
both of these research directions. 

Many models predict a diffuse flux of high energy neutrinos from active galactic nuclei 
and other astrophysical sources. The "Astrophysics Underground" portion of this thesis 
describes a search for this neutrino flux performed by looking for extremely high energy 
upward-going muons using the Super-Kamiokande detector, and comparing the observed 
flux to the expected background. We use our results to to set an upper limit on the diffuse 
neutrino flux from astrophysical sources. 

In addition to using particles to do astronomy, we can also use the universe itself as 
a particle physics lab. Cosmology provides new insights that could never be observed in 
terrestrial laboratories. The "Particle Physics in the Sky" portion of this thesis focuses on 
extracting cosmological information from galaxy surveys. 

To overcome technical challenges faced by the latest galaxy surveys, we produced a 
comprehensive upgrade to mangle, a software package that processes the angular masks 
defining the survey area on the sky. We added dramatically faster algorithms and new 
useful features to this software that are necessary for managing complex masks of the Sloan 
Digital Sky Survey (SDSS) and will be invaluable for future surveys as well. 

With this software in hand, we utilized galaxy clustering data from SDSS to investigate 
the relation between galaxies and dark matter by studying relative bias, i.e., the relation 
between different types of galaxies. If all galaxies were perfect tracers of dark matter, 
different subpopulations would trace each other perfectly as well. However, separating 
galaxies by their luminosities and colors reveals a complicated picture: red galaxies are 
clustered more strongly than blue galaxies, with both the brightest and the faintest red 
galaxies showing the strongest clustering. Furthermore, red and blue galaxies tend to occupy 
different regions of space, effectively introducing an element of stochasticity (randomness) 
when modeling their relative distributions. In order to make precise measurements from 
the next generation of galaxy surveys, it will be essential to account for this complexity. 

Thesis Supervisor: Max Erik Tegmark 
Title: Associate Professor of Physics 



3 



4 



Acknowledgments 



I would like to extend my heartfelt thanks: 

To my first teachers, my parents Barb and Rick Swanson. To my dad for everything 

from counting Cheerios to help with calculus homework, and for teaching me that physicists 
answer questions in the language of math, even when their 8-year-old daughters can't bear 
to see another equation scribbled on a napkin, and to my mom for her endlessly enthusiastic 
support of everything I've ever done. 

To my most recent teachers, my advisors Max Tegmark and Kate Scholberg, for being 
such great role models, for always sticking up for me, and for showing me two such vastly 
different yet equally inspiring ways to be a scientist. 

To all the great teachers I've had in between. 

To my late grandmother Helen Hage - it was worth it, even though I'm not a boy. 

To all of the Dr. Swanson physicists before me - my dad, my uncle Charlie Swanson, 
and my grandfather Bob Swanson, thanks to whom I always thought getting a PhD in 
physics was a common thing to do. 

To the rest of the Swanson(/01son/McLaughlin) clan, for doing other things so that we 
are not all huge nerds. 

To the Crosby/Newman/ Allen/etc. clan, especially to my grandfather-in-law Da (known 
to the rest of the world as Bob Newman) for taking a genuine interest in what I do and for 
my favorite definition of a neutrino ever. 

To the Techcr-Boston crew and various other dear friends in the area, including Wedge 
Cheung, Robin Friedman, Mike Fucrstman, Todd Schuman, Dave Guskin, Mick Garvey, 
Marlena Fecho, Brian and Sarah Bairstow, Craig Chu, Charley Mills, Pat Codd, Dave and 
Jess Tytell, Kate Jensen, and Dan Recht. The Thursday-night-dinners-on-Tuesdays-on- 
Wednesdays-on-Thursdays are among my fondest memories of my time in Boston. 

To all my fellow astrograds, and everyone else I've had the pleasure of knowing at MIT, 
especially Miriam Krauss, Jake Hartman, Judd and Cassie Bowman, Ed Boyce, Robyn 
Sanderson, Ben Cain, Mike Stevens, Yi Mao, Bonna Newman, Julie Millane, Jessie Shelton, 
Tom and Liz Pasquini, Will Fox, and Sarah-Jane White. Having such awesome people with 
whom to share the pain made grad school infinitely more fun. 

To my fellow 37-602 denizens, most especially Josh Carter and Ryan Lang, who were 
in it with me for the long haul, but also to Judd Bowman, Sam Conner, Chris Williams, 
Leslie Rogers, and (albeit briefly) Sarah Vigeland. It wouldn't have been the same without 
dodging Nerf footballs, green tea, Penzey's spice rub, indoor paper baseball, cutthroat 
games of Facebook Scrabble and World Traveler, multiwing bats, SoftLips pranks, paper 
airplane contests, and countless other diversions, not to mention the constant support. 

To Gideon Koekoek, virtual officemate extraordinaire, who, thanks to the magic of teh 
intarnets, was just as distracting and just as inspiring as his real-life counterparts all the 
way from Amsterdam. Words can't do it justice, so here are some *doctoral hugs* instead. 

To all of Walking Fish and Blacker Hovse, for being the most amazing friends a nerdy 
girl could ask for, especially to Katie Romportl, my first partner-in-crime, to Mike Borchert 
and Matt Kuzma for throwing pennies at my window late at night in order to sneak into my 
room to chat about quantum mechanics, to Mike Jorgenson for fractals and ensanguination, 
to Joe Carroll for frisbee and folk songs, to Nate Austin for castles and zepplins, and to 
Katharina Kohler for always understanding exactly what I'm going through. 

And most of all, to my husband Tim Crosby, for being the reason I come home every 
night, and for being the home I'll always come home to. Swarling swar swar, swarling swar. 



5 



The work presented in chapters O O [6] and ^4.31 was ah done in collaborations - here I 
would like to acknowledge my co-authors and detail my own role in each of these projects. 

For chapter [3l I did the analysis and write-up primarily by myself, with valuable guidance 
from Kate Scholberg, Alec Habig, Shantanu Desai, and Jodi Cooley. The other members 
of the Super-Kamiokande collaboration helped via building the detector, keeping it run- 
ning, and providing insightful comments on drafts of the paper. I also gratefully gratefully 
acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super- 
Kamiokande experiment has been built and operated from funding by the Japanese Ministry 
of Education, Culture, Sports, Science and Technology, the United States Department of 
Energy, and the US National Science Foundation. 

For chapter El I added the new pixelization algorithms described here to the mangle 
software which was originally written by Andrew Hamilton, and did the bulk of the write-up 
as well. The HEALPix algorithms were done primarily by UROP student Colin Hill. Colin 
and I also worked together to create the new mangle website - Colin deserves a special 
thank you for pulling an all-nighter with me the night before the website went live. I also 
thank Michael Blanton for providing the SDSS DR6 VAGC mask, Krzysztof Gorski and 
collaborators for creating the HEALPix package, and the WMAP team for making their 
data public via the Legacy Archive for Microwave Background Data Analysis (LAMBDA; 
http : / /lambda . gsf c . nasa . gov) . . Support for LAMBDA is provided by the NASA Office 
of Space Science. 

For §4.31 which summarizes Tegmark et al. ( 20061 ) . I contributed primarily by using my 
new version of MANGLE to create angular masks for each of the angular subsamples used in 
the analysis. I also helped with proofreading and editing the original paper. 

For chapter [U I did the analysis and write-up, with guidance from Max Tegmark. Idit 
Zehavi provided careful editing and insightful suggestions, and Michael Blanton provided 
the SDSS VAGC-DR5 data. I also wish to thank Daniel Eisenstein, David Hogg, Taka Mat- 
subara, Ryan Scranton, Ramin Skibba, and Simon White for helpful comments. Funding 
for the SDSS has been provided by the Alfred P. Sloan Foundation, the Participating In- 
stitutions, the National Science Foundation, the U.S. Department of Energy, the National 
Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck 
Society, and the Higher Education Funding Council for England. The SDSS Web Site is 
[http : //www, sdss . org , The SDSS is managed by the Astrophysical Research Consortium 
for the Participating Institutions. The Participating Institutions are the American Museum 
of Natural History, Astrophysical Institute Potsdam, University of Basel, Cambridge Univer- 
sity, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, 
the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, 
the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics 
and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), 
Los Alamos National Laboratory, the Max-Planck- Institute for Astronomy (MPIA), the 
Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State 
University, University of Pittsburgh, University of Portsmouth, Princeton University, the 
United States Naval Observatory, and the University of Washington. 

This work was supported by a National Defense Science and Engineering Graduate Fel- 
lowship, a Bruno Rossi Fellowship, NASA grant NNG06GC55G, NSF grants AST-0134999, 
0607597 and 0708534, the Kavli Foundation, and fellowships from the David and Lucile 
Packard Foundation and the Research Corporation. 



6 



Contents 



List of Figures 
List of Tables 

1 Connecting particle physics and astrophysics 

1 Astrophysics Underground 

2 Neutrino astronomy 

2.1 History of multi-messenger astronomy 

2.2 Neutrinos as astronomical messengers 

2.3 Astrophysical production mechanisms 

2.3.1 Diffusive shock acceleration 

2.3.2 Active galactic nuclei 

2.3.3 Other sources 

2.4 Neutrino telescopes 

2.4.1 How can we detect the astrophysical neutrino flux? 

2.4.2 Super-K 3jS Si neutrino telescope 

3 Search for diffuse astrophysical neutrino flux using ultra-high energ 
upward-going muons in Super-Kamiokande I 

3.1 Introduction 

3.2 The Super-Kamiokande Detector 

3.3 Event Selection 

3.3.1 Outer Detector Linear Fit 

3.3.2 Selection Cuts 

3.4 High-Energy Isotropic Monte Carlo 

3.4.1 Resolution and Efficiency of Event Selection 

3.4.2 Ultra-High-Energy Fraction 

3.5 Flux Calculation 

3.6 Expected Atmospheric Background from Monte Carlo 

3.7 Analytical Estimate of Expected Muon Flux 

3.7.1 Method for Calculating Muon Flux 

3.7.2 Analytical Estimate of Expected Background 

3.8 Upper Limit for Muon Flux from Cosmic Neutrinos 

3.9 Conclusions 



7 



II Particle Physics in the Sky 



55 



Doing cosmology with galaxies 

4.1 Cosmology Basics 

4.1.1 What is cosmology? 

4.1.2 The standard cosmological model 

4.1.3 Zeroth order cosmology 

4.1.4 First order cosmology 

4.1.5 Higher order cosmology 

4.1.6 What does this teach us about particle physics? 

4.2 Tools of cosmology 

4.2.1 Combining cosmological probes 

4.2.2 Galaxy surveys 

4.2.3 Cosmic microwave background 

4.2.4 Other probes 

4.3 Example: combining WMAP and SDSS 

4.3.1 Power spectrum estimation 

4.3.2 Cosmological parameter estimation 

4.4 Issues with using galaxies as cosmic tracers 

4.4.1 Managing angular masks of galaxy surveys . . 

4.4.2 Relative galaxy bias 



Methods for rapidly processing angular masks of next-generation galaxy 
surveys | 

5.1 Introduction 

5.2 Mangle terminology 

5.3 Speedup: pixelization 

5.3.1 Pixelization concept 

5.3.2 Pixelization schemes 

5.3.3 Pixelization algorithm 

5.3.4 Speed trials 

5.4 Unification with HEALPix and other pixelized tools [iol 

5.4.1 Importing HEALPix maps into mangle 

5.4.2 Exporting polygon files as HEALPix maps 

5.5 Summary 



104 
10£ 



6 SDSS galaxy clustering: luminosity & color dependence and st o chast icit v lOol 



6.1 Introduction [io^ 

6.2 SDSS Galaxy Data 

6.3 Analysis Methods 

6.3.1 Overlapping Volume-Limited Samples 

6.3.2 Counts-in-Cells Methodology 

6.3.3 Relative Bias Framework 

6.3.4 The Null-buster Test 

6.3.5 Maximum Likelihood Method 

6.4 Results 

6.4.1 Null-buster Results 

6.4.2 Likelihood Results 



8 



6.5 Conclusions 

6.5.1 Biasing results 

6.5.2 Implications for galaxy formation 

6. A Consistency Checks 

6.A.1 Alternate null-buster analyses 

6.B Uncertainty Calculations 

6.B.1 Jackknife uncertainties for null-buster analysis 
6.B.2 Likelihood uncertainties 



7 Conclusions 
Bibliography 



9 



10 



List of Figures 



1- 1 Flow diagram of connections between particle physics and astrophysics. . . 

2- 1 Geometry of diffusive shock acceleration 

3- 1 Schematic drawing of the Super-K detector 

3-2 OD-based muon trajectory fit applied to example MC muon events 

3-3 Resolution of OD-based fit on determining track parameters for events from 

the high energy isotropic MC 

3-4 Angular resolution of OD-based fit for events from the high energy isotropic 
MC 

3-5 Distribution of the time difference between OD entry and ID entry clusters 
for UHE downgoing muon events from data, compared to the distribution for 
the UHE events in the high-energy isotropic MC 

3-6 Efficiency of our data reduction on upward, throughgoing, > 1.75 x 10^ pe 
muons 

3-7 Fraction of events k {Efj_) that deposit > 1.75 x 10^ pe in the ID for the 100 
yr atmospheric MC, the high-energy isotropic MC, and a power-law fit to the 
results from the latter 

3-8 The number of upward-going muon events above energy threshold in- 
duced by atmospheric neutrinos, from our analytical calculatoin and from 
the 100 year atmospheric MC 

3-9 Upper limits from this analysis on muon (/x"'' + flux above energy thresh- 
old -E™™, compared to various model fluxes 

3- 10 Approximate upper limits from SK-I on astrophysical neutrinos {v^ + P^). . 

4- 1 The curvature of space illustrated for two-dimensional universes 

4-2 A slice through the CfA redshift survey 

4-3 A slice through 2dFGRS 

4-4 A slice through SDSS 

4-5 Temperature anisotropics in the CMB from WMAP5 

4-6 CMB angular power spectrum from WMAP5 

4-7 The angular distribution of the SDSS DR4 luminous red galaxies showing 

the seven angular subsamples analyzed 

4-8 Measured power spectra for the full LRG and main galaxy samples 

4- 9 95% constraints in the {ujd, fu) plane 

5- 1 Speed trials for a series of portions of the SDSS DR5 mask with and without 

pixelization 

5-2 A portion of the SDSS DR5 angular mask 



11 



94 
97 
97 



5-3 A cartoon illustrating the process of balkanization (a) with no pixelization, 
(b) with pixelization 

5-4 The full sky pixelized with the simple pixelization scheme 

5-5 The full sky pixelized with the SDSSPix pixelization scheme 

5-6 Left: portion of SDSS mask from Fig. 15-21 pixelized to a fixed resolution with 
the simple pixelization scheme. Right: the same mask pixelized with the 
simple pixelization scheme using the adaptive resolution method lOOl 

5-7 Time required to process the full SDSS DR5 mask with different choices of 



pixelization schemes and methods llOl 

5-8 Left: Portion of SDSS mask as shown in Fig. 15-21 Right: Portion of SDSS 



mask from Fig. 15-21 as approximated by HEALPix pixels 1051 



5-9 Top: The SDSS DR5 completeness mask, rasterized and plotted using HEALPix 
routines. Middle: The final 2dFGRS completeness mask, rasterized and plot- 
ted using HEALPix routines. Bottom: CMB temperature difference map 
measured by WMAP channel 4 lod 



6-1 Histogram of the comoving number density of the volume-limited samples. . 1131 
6-2 Color-magnitude diagram showing the number density distribution of the 

galaxies in the volume- limited samples 1151 

6-3 Galaxy distributions plotted in comoving spatial coordinates for a radial slice. 



118 



118 



showing the four different types of pairwise comparisons made lid 

6-4 The SDSS DR5 angular mask pixelized at the four different resolutions used 

to partition the survey into cells 

6-5 A radial slice of the SDSS survey volume divided into cells 

6-6 Null-buster results for pairwise comparisons 

6-7 Relative bias 6rei between pairwise samples 

6-8 Luminosity dependence of bias for all, red, and blue galaxies 

6-9 Comparison to previous results for the luminosity dependence of bias for all. 



128 



129 



131 



red, and blue galaxies 1321 

6-10 The best-fitting values of the relative cross-correlation coefficient rj-^i between 



136 



pairwise samples Il35 

6-11 Comparison of relative cross-correlation coefficient rj-^i between red and blue 

galaxies as measured with different techniques 

6-12 Null-buster results for randomly split samples 

6-13 Comparison of jackknife and generalized uncertainties on 6rei from the 

null-buster analysis 

6-14 Typical contour plots of A (21n£) for volume V5 for three different resolutions 



13S 



142 



143 



12 



List of Tables 



3.1 Visual scan of candidate upward-going muon events 

3.2 Fraction of high-pe events k (E^) in high energy isotropic MC 

3.3 The flux of ultra-high energy upward-going muons as observed by SK-I. . . 

3.4 Systematic uncertainties in atmospheric neutrino background 

3.5 Confidence intervals for the upward-going muon flux due to neutrinos from 
astrophysical sources 

3.6 Approximate upper limits from SK-I on astrophysical neutrinos (f^ -|- P^). . 



6.1 Summary of cuts used to create luminosity-binned volume- limited samples. 

6.2 Overlapping volumes in which neighboring luminosity bins are compared. . 

6.3 Number of galaxies in each sample being compared 



36 
41 
42 
44 



49 
51 



4.1 Cosmological parameters measured from WMAP and SDSS LRG data. ... \8l 
5.1 Definitions of Terms, in Alphabetical Order [o^ 



114 



114 



117 



13 



14 



Chapter 1 



Connecting particle physics and 
astrophysics 

Meg looked about. Ahead of her was a tremendous rhythmic swirl of wind 
and flame, but it was wind and flame quite different from the cherubim's; this 
was a dance, a dance ordered and graceful, and yet giving an impression of com- 
plete and utter freedom, of ineffable joy. As the dance progressed, the movement 
accelerated, and the pattern became clearer, closer, wind and fire moving to- 
gether, and there was joy, and song, melody soaring, gathering together as wind 
and fire united. 

And then wind, flame, dance, song, cohered in a great swirling, leaping, 
dancing, single sphere. 

Meg heard Mr. Jenkins's incredulous, "What was that?" 
Blajeny replied, "The birth of a star." 

Mr. Jenkins protested, "But it's so small I can hold it in the palm of my 
hand." And then an indignant snort, "How big am I?" 

"You must stop thinking about size, you know. It is both relative and 
irrelevant." 

— Madeleine L'Engle, A Wind in the Door 

One of the most remarkable things about physics is that it can be applied over such 
a vast range of size scales. Some physicists focus on understanding our world on almost 
inconceivably tiny scales, studying theories of fundamental building blocks of matter as 
small as 10~^^ meters, a factor of 10~^° smaller than the diameter of a proton. Other 
physicists study our world at equally inconceivably huge scales, looking at the most distant 
parts of our universe that it is possible to sec, over 10^^ meters away. These largest and 
smallest scales in physics are separated by over 60 orders of magnitude. 

Perhaps even more astounding is the fact that these two extreme regimes of physics are 
deeply connected: learning about how our universe behaves on the very largest scales can 
teach us things about elementary particles that could never be discovered in a terrestrial 
laboratory, and studying particles coming from space can also teach us about properties of 
distant astronomical objects. Over the past several decades, physicists have gained many 
great insights about both the macroscopic and microscopic regimes by taking advantage of 
this connection. The link between particle physics and astrophysics, which is the focus of 
this thesis, is but one example of a fruitful interaction between traditionally separate fields 
of science - there is a tremendous amount of exciting work being done at interfaces between 
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Figure 1-1 Flow diagram of connections between particle physics and astrophysics. 



other fields as well. Some examples include applying statistical physics to biological systems, 
quantum computation research that aims to use quantum physics principles to revolutionize 
computer science, and neuroscience research integrating biology, psychology, and computer 
science. 

Doing science at these interfaces is no easy task - integrating expertise across differ- 
ent fields generally leads to certain challenges. Scientists on either side of the fence will 
have different academic cultures, traditions, and methodologies, and they frequently do not 
mix well. For example, concerns have been raised that the trend towards using astronomy 
to st udy fundamen tal physics could ultimately be damaging to the astronomical commu- 
nity: IWhitd (|2007l ) notes that "By uncritically adopting the values of an alien system [of 
high-energy physics traditions] , astronomers risk undermining the foundations of their own 
current success and endangering the future vitality of their field." However, the scientific 
potential of the crossover is large enough that such obstacles are well worth overcoming. As 
Kolbl (l2007l ) says, the ties between astronomy and particle physics "can either encourage 



or choke creativity ... It is up to us to choose wisely." If these two scientific communi- 
ties intelligently navigate the cultural interface, the science of both fields will be greatly 
enriched. 

The nature of the particle-physics-astrophysics connection used in this thesis is illus- 
trated in Fig. 11-11 The traditional means of doing astrophysics is to look at the light emitted 
from distant astronomical objects with a telescope and use this to infer the physics of these 
objects. On the other hand, traditional particle physics is done by measuring the properties 
of particles passing through a detector and using this to infer the nature of the particles 
themselves or the interactions that produced them. 

One means of connecting astrophysics and particle physics is to use particles as astro- 
nomical messengers, i.e. to look for particles other than photons coming from astronomical 
objects to get a more complete picture of the physical system. We dub this "Astrophysics 
Underground" because when the particle messenger is a neutrino, the "telescopes" used for 
this endeavor are typically large particle detectors located deep underground to shield from 
background radiation. Essentially this is using the tools of particle physics to do astronomy. 

We can also turn this connection in the other direction and use the tools of astronomy to 
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do particle physics, which we dub "Particle Physics in the Sky." One particularly exciting 
area in which this can be applied is cosmology, the study of our universe as a whole. By 
doing astronomical observations cosmologists have built up a consistent theory of how our 
universe has evolved that turns out to depend on the properties of elementary particles and 
has also provided evidence for physics beyond the standard model of particle physics. Thus 
by studying the our universe at the largest scales we can also learn about the smallest. 

The rest of this thesis is organized as follows: an exploration into doing "Astrophysics 
Underground" using neutrinos is given in chapters [2] and [3l Chapter [2] provides an overview 
of the field of neutrino astronomy, and chapter [3] describes a search for high energy neutrinos 
coming from astrophysical sources using the Super-Kamiokande neutrino detector. The 
following three chapters turn the connection arrow in the other direction and focus on 
"Particle Physics in the Sky" by doing cosmology with galaxy surveys. Chapter U] gives a 
general introduction to the field of cosmology with a particular focus on the tool of galaxy 
surveys and what this can tell us about particle physics. In chapter [5l we describe a set of 
computational techniques for managing the complex data of modern galaxy surveys, and 
in chapter [6] we use data from the Sloan Digital Sky Survey to study how the clustering of 
galaxies depends on and their luminosities and colors and the implications of this for using 
galaxies as cosmological tools. Finally, we conclude in chapter [7| with a summary of what 
we have learned, some ideas for future research directions, and reflections on the nature of 
doing science at this fascinating interface between astrophysics and particle physics. 
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Part I 

Astrophysics Underground 
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Chapter 2 

Neutrino astronomy 



2.1 History of multi- messenger astronomy 

Historically, the study of astronomy has been the study of photons coming from the sky. 
Early astronomers learned a tremendous amount observing only the visible- wavelength pho- 
tons that could be seen with the naked eye. Since then, revolutions in astronomy have come 
about through expanding our ability to detect different types of astronomical messengers. 

The invention of the telescope in the early seventeenth century enabled the collection of 
much larger numbers of photons, opening a window on fainter sources. The twentieth cen- 
tury saw an explosion of astronomical information as techniques were developed to observe 
photons all across the electromagnetic spectrum from radio waves to gamma rays, leading 
to the discovery of remarkable astrophysical phenomena such as pulsars and gamma ray 
bursts. 

Detecting particles other than photons from astrophysical sources provides yet another 
window on the universe. The discovery of cosmic rays in 1912 was the first step in this 
endeavor. Cosmic rays are charged particles - protons, electrons, and atomic nuclei - origi- 
nating from cosmic sources. Studies of cosmic rays have led to insights about astrophysical 
processes such as the evolution of galactic magnetic fields and particle acceleration within 
remnants of supernova explosions. High-energy cosmic rays penetrating the atmosphere 
also provide a natural particle collider - this was an instrumental tool in particle physics 
as well: for example, positrons, muons, and pions were all discovered in cosmic ray studies. 
Today the astrophysical origin and acceleration mechanisms of the highest energy cosmic 
rays is one of the primary open questions in particle astronomy. 

2.2 Neutrinos as astronomical messengers 

Another exciting step towards using particles as astronomical messengers is to look for 
neutrinos coming from cosmic sources. In terms of providing astrophysical information, 

neutrinos have some advantages over photons and charged particles. Firstly, they are not 
readily absorbed by intervening matter, so they can act as probes of extremely high-density 
environments such as the central cores of stars, active galactic nuclei, or supernovae from 
which photons could never escape. Secondly, they are not deflected by cosmic magnetic 
fields, so unlike charged particles, they always point directly back to their source. 

The advantage of low absorption, however, is the biggest disadvantage as well - they 
tend to pass through our detectors unseen just as readily as they pass through dense ma- 
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terial at the source. Thus neutrino astronomy demands sophisticated detection techniques, 
and as such, this field is still in its inf ancy. So far, neutriiios ha ve only been detected 
from two astrophysical sources: the su n (jPavisI I2OO3I : iKoshibal liooi l and supernova 1987a 
(jHirata et al.l 119871 : iBionta et al.l 119871 ). However, the insights gained into both astronomy 
and particle physics from just these two neutrino sources have been tremendous. 

The work of Raymond Davis and others on detecting neutrinos from the sun gener- 
ated the famous "solar neutrino problem", whereby the observed flux of solar neutrinos 
was only about a thi rd of what w as expected based on astrophysical models of energy pro- 
duction in the sun (|Dayislll99d l. This called both the experimental measurements and 
the standard solar models into question, and remained a controversy for decades. Its ul- 
timate resolution actually came from a revolution in particle physics: measurements from 
the Super-Kamiokande neutrino detector (Super-K; Fukuda et al. 19981 ) provided conclusive 
evidence that neutrinos oscillate from one type to another and therefore must have mass, 
in conflict with the standard model of particle physics which treats them as massless. This 
explained the solar neutrino problem via some fraction of electron neutrinos produced in 
the sun oscillating to muon or tau neutrinos by the time they reach the Earth and thus 
going undetected. Ultimately, observations of solar neutrinos have both solidified our un- 
derstanding of astrophysical processes in the sun and provided the first concrete evidence 
for physics beyond the standard model. 

The 20 neutrinos detected from supernov a 1987a by the Kamiokande-II and 1MB de- 
tectors ( Hirata et al. 1987 : Bionta et al. 19871 ) similarly provided insights into both particle 
physics and astrophysics. The flux observed is consistent with astrophysical models of 
core-collapse supernovae in which nearly all of energy from the collapse is emitted in neu- 
trinos. Additionally, the observation that neutrinos of different energies all arrived nearly 
simultaneously provides a clue to the nature of the neutrinos themselves: it means that 
the neutrinos of different energies all travel at nearly the speed of light, which implies 
that the electron neutrino rest mass must be quite small. Thus the neutrinos observed 
from the supernova can be used to set an upper limit on the mass of the electron neutrino 
( Kernan and Krausi 19951 ). 



The next frontier of neutrino astronomy is to detect astrophysical neutrinos at high 
energies (GeV-PeV) - several experimental efforts to do this are currently underway. Such 
high-energy neutrinos will open yet another new window on the universe - they will shed 
light on astrophysical processes in extreme env ironments and will also p rovide clues to the 
origins of the highest energy cosmic rays. (See Halzen and Hooper 20021 for more details.) 



2.3 Astrophysical production mechanisms 
2.3.1 Diffusive shock acceleration 

The general mechanism for producing high-energy neutrinos is the following: first, protons 
or other nuclei are accelerated by shocks in an astrophysical plasma. Then the accelerated 
hadrons interact with photons in a surrounding radiation field and produce charged pions 
that decay into high-energy neutrinos: 
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These interactions can also produce neutral pions, which generate high-energy gamma 
rays: 
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The proton acceleration mechanism at work h ere is known as diffusive shock acceleration, 
also known as the first-order Fermi mechanism ( Fermi 19491 ) . The difference in velocities 
on either side of a shock in a magnetized plasma can transfer kinetic energy to particles 
as they cross the shock multiple times. The basic geometry of a plane shock is illustrated 
in Fig. 12-11 This mechanism generically predicts a power-law spectrum of particle energies 
with dNidE oc E~°' where a ~ 2. We give a rough deriv a tion o f this below following 
Gaisseil (|l99nl ): for a more detailed treatment see ISchlickeiseil jiooi). 

The basic equations governing the behavior of shocks can be derived from conserva- 
tion of mass, momentum, and energy across the shock boundary. For a non-relativistic 
hydrodynamic plane shock, these conditions are given by 



0, (2.4) 
P1-P2, (2.5) 

0, (2.6) 

where 1 subscripts denote conditions upstream of the shock and 2 subscripts denote down- 
stream conditions, u is the component of the velocity normal to the shock front in the rest 
frame of the shock, p is the density, P is the pressure, and w is the specific enthalpy, given 
by u; = CpT for temperature T and specific heat at constant pressure Cp. For an ideal gas, 
w = 7P/ ((7 — 1) p) where 7 is the ratio of specific heats at constant pressure and constant 
volume; 7 = 5/3 for a monatomic ideal gas. In a magnetized plasma the magnetic fields 
will also contribute to these conditions - here we assume either that the magnetic field is 
parallel to the normal of the shock plane or that the magnetic fields are small enough to be 
negligible for the continuity conditions. 

In order for a shock to form, ui must be greater than the sound speed Cgi = \/ 7-P1 / pi ■ 
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Figure 2-1 Geometry of diffusive shock acceleration, shown in the rest frame of the material 
upstream of the shock. The trajectory of an accelerated particle is shown in blue/grey. 
Scattering centers (i.e. magnetic field irregularities) are represented as circles. 



The Mach number M of the shock is defined as M = ui/csi, and by using equations (j2.4p . 
(|2.5p . and (12.61) . and much clever algebra, the ratio of the upstream to downstream velocities 
is given by 

«i ^ (7 + 1) 

U2 M2(7-l) + 2- ^ > 

The energy of a charged particle will increase as it diffuses across the shock from the 
upstream medium to the downstream medium and back, as illustrated by the trajectory in 
Fig. 12-11 To calculate the energy change, we start with a relativistic {E ~ pc) population of 
particles with velocities that are isotropic in the frame of the upstream medium. A particle 
with energy Ei that crosses the shock at an angle 9i in the upstream frame will have an 
energy 

E[ = T El {I- 13 cos Oi) (2.8) 

in the downstream frame, where [5 = (ui — U2) /c and F is the corresponding Lorentz factor. 
Primes denote quantities in the downstream frame; unprimed quantities are in the upstream 
frame. Averaging over an isotropic population gives (cos^i) = —2/3 for particles crossing 
the shock. 

In the downstream medium, the particle undergoes elastic scattering off of irregularities 
in the magnetic field generated by low-frequency magnetohydrodynamic waves. As the 
particle diffuses, its average motion will coincide with that of the downstream medium, so 
we now have an isotropic population in the downstream frame. To a first approximation, 
these scatterings do not cause the particle to lose energy, so a particle that is scattered back 
into the upstream medium at an angle 62 will have energy E2 = E[. Transforming back 
into the upstream frame gives 

^2 = rS^l + /?cos02) , (2.9) 
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and since the particles have isotropized relative to the downstream medium, (cos02) can 
be obtained by averaging over an isotropic population as well, yielding (cos ^2) — 2/3. 
Plugging equation (j2.8p in for E'2 in equation (j2.9p and averaging over the angular factors 
yields the following formula for the fractional energy change in one "encounter" (diffusing 
from upstream to downstream and back): 



1 + 1/3 + 1/3^ 

l-/?2 

« (2.10) 

where the last approximation holds for a non-relativistic shock with /? <C 1. 

Once in the downstream medium, the particles have some probability to escape rather 
than get scattered back into the upstream medium. We can estimate this probability by 
considering the flux of particles through an imaginary plane one scattering length away from 
the shock in the downstream medium. Particles that pass through this plane back in the 
upstream direction will get swept up by the shock, and particles that pass through towards 
the downstream direction can escape. Working in the frame of the downstream medium, 
particles will pass back through the shock if ccos9 < —U2- Thus the flux upstream is 

(ccos 6 + U2) dcos 9 

_ (C - U2f 

2c 

and the flux downstream is 
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The probability for a particle to escape is then 



/(ccos 9 + U2) dcos 9 
-uojc 
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(2.11) 



for U2 <^ c. 



The energy spectrum produced in this scenario - where a particle gains fractional energy 
e and has a probability Pesc of escaping with each encounter ~ can be calculated as follows: 
if particles initially have energy £"0, the energy after m encounters is En = Eo{l + e)^, so 
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the number of encounters n needed to reach energy E is 



Tl — 

log(l + e)' 

The probabihty of remaining after m encounters is (1 — Pcsc)™', so the number of particles 
with energy greater than E is thus 

oo 
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a - 1 = , (2.12) 
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The differential energy flux is thus dN/dE oc Plugging in the results of equa- 

tions (|2.1U|) and (j2.11|) into equation (|2.12|) and using the approximations that ui and U2 <C 
c, 

3 



a « 1 + 



U2 



Now we can insert equation (j2.7p for ui/u2^ use 7 = 5/3 for a monatomic ideal gas, and 
take the limit of a strong shock where M ^ 1, yielding 

a«2 + A. (2.13) 

Thus, dN/dE ~ i?"^ for first order Fermi acceleration. 

Calculating the details of the high-energy neutrino spectrum of a particular type of astro- 
physical source requires a more sophisicated treatment that relaxes some of the assumptions 
we have made here and accounts for appropriate energy loss mechanisms as well as how 
the pion energy is distributed among the decay products in a particular context. However, 
since most astrophysical models invoke this first-order Fermi mechanism, the prediction of 
a spectrum roughly proportional to is a generic feature. As such, d^^y/dEy oc E~'^ is 
frequently used as a fiducial model for an astrophysical neutrino flux - we make use of this 
convention in chapter [3l 



2.3.2 Active galactic nuclei 

Astrophysical systems that produce high energy neutrinos via the mechanism outlined in the 
previous section must have two features: a shocked, magnetized plasma that can accelerate 
protons to very high energies, and a sufficiently intense radiation field to allow for photo- 
pion production. These conditions are expected to be present in a number of astrophysical 
systems, including gamma ray bursts and active galactic nuclei (AGNs). The expected 
signal from AGNs is the strongest for the energy range to which Super-K is sensitive, so 
they are the primary astrophysics target for the study discussed in chapter [3l 
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AGNs are supermassive black holes accreting matter at the centers of galaxies. They are 
the most luminous persistent sources in the universe and emit radiation over a vast range 
of the electromagnetic spectrum, from radio waves to gamma rays. There are two locations 
where an AGN could produce neutrinos: in the core or in the jets. We give a rough outline 
of these processes here; further details can be found in lGaisser et al. I (1 19951 ). 



In models of neutrino production within AGN cores, a shock is formed in an approxi- 
mately spherical shell around the black hole, where the ram pressure of the infalling material 
is balanced by the radiation pressure of the emission. The region inside the shock is opti- 
cally thick and provides ample radiation density for p7 collisions. This process generates 
high-energy neutrinos from quite close to the black hole, thus probing a region from which 
photons cannot escape. A detailed model of neutrino production in AGN cores can be 
found in Stecker and Salamon (199^), along with a prediction for the diffuse flux of neutri- 



nos expected from integrating the contributions from all AGN cores over the history of the 
universe. 

Neutrinos can also be produced within AGN jets. Jets are relativistic outflows of charged 
particles collimated along the black hole's rotation axis. The mechanism that powers the 
jets is not yet well-understood in detail, but the basic idea is that a varying magnetic field 
linked to the rotation of the black hole or the accretion disk generates an electric field 
that can accelerate charged particles. In the case of jets, the shocks needed for first-order 
Fermi acceleration to extremely high energies are formed by inhomogeneities in the matter 
traveling along the jet. The required photon field is provided by synchrotron radiation from 
electrons spiraling around magnetic field lines. A model of the predicted diffuse neutrino 



flux from AGN jets can be found in lMannheim et al.l (|200lh . 



It is not known whether jets consist mainly of electrons and positrons or if there is a 
hadronic component as well. The observed gamma rays coming from jets could be explained 
by either leptonic or hadronic processes. In the leptonic process, gamma rays are produced 
via low-energy synchrotron photons gaining energy via inverse Compton scattering with 
the accelerated electrons - this is known as synchrotron-self-Compton radiation. Alterna- 
tively, accelerated hadrons could produce gamma rays via the decay of neutral pions, as in 
equation (j2.3p . Since neutrinos are only produced in hadronic processes, observation of a 
neutrino flux from these sources would help distinguish between these models. The question 
of whether hadrons are accelerated to extremely high energies by AGNs is also an important 
question for determining the origin of ultra- hi gh energy cosmic rays. Recent results from 
the Auger Observatory ( Abraham et al. 200?! ) indicate that the arrival directions of the 



highest energy cosmic rays are statistically correlated with the locations of nearby AGNs 
- this is a tantalizing hint that AGNs could be the answer to the long-standing puzzle of 
ultra-high energy cosmic rays. A definitive obser vation of a neutrino flux associ ated with 
AGNs would lend strong support to this picture ( Halzen and O'Murchadha 20081 ). 



2.3.3 Other sources 

In addition to AGNs, there are many other possible astrophysical sources of neutrinos. 
For example, gamma ray bursts (GRBs) - the highest-energy explosions in the universe 
- are likely to provide t he conditions needed for dif fusive shock acceleration and photo- 
pion production as well ( Waxman and Bahcall 19991 ). There are also a number of more 



exotic phenomena that could produce high-energy neutrinos, including decays of topolog- 
ical defects such as cosmic strings or magn etic monopol es and self-annihilation of weakly 
interacting dark matter particles (WIMPs) (jStastolliood ). 
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Another source of high-energy neutrinos is hnked to the propagation of ultra-high en- 
ergy cosmic rays through the universe: protons with energies above ~10^° eV will un- 
dergo photo-pion production as in equation (|2.3p via collisions with cosmic microwave 
background (CMB) photons (see ^4.2.3I for more abo ut the CMB). This is known as the 
GZK cutoff (Greisen 1966 : Zatsepin and Kuzmin 19661 ) . and recent cosmic ray experiments 
(jSokolsky and Thomson! 120071 ) claim to have observed the predicted break in the cosmic 
ray spectrum due to this effect. Thus there should be a correspond ing neutrino flux from 



the decay of the pions produced - experiments such as ANITA f^Barwick et al 
attempting to detect these neutrinos. 



20061 ) 



are 



2.4 Neutrino telescopes 



2.4.1 How can we detect the astrophysical neutrino flux? 

The flux of high-energy neutrinos from active galactic nuclei is quite small and the proba- 
bility of a neutrino interacting with material in a detector is also quite low. These two basic 
facts point to the need for an extremely large collecting area. Modern neutrino "telescopes" 
are thus actually large naturally-occurring volumes of ice or water instrumented with pho- 
tomultiplier tubes (PMTs). Neutrinos that interact in or near the detector volume produce 
charged particles that emit Cerenkov light as they pass through the water or ice; this light 
is then recorded by the PMTs. 

In order to make a significant observation of astrophysical neutrin o flux, it is expected 
that we will need a detector with a volume on the order of Ikm^ ( Halzen and Hooper 



2OO2I ). There are two main efforts underway aiming towards this goal: IceCube, a detector 
currently under construction at the South Pole which uses the Antarctic ice shelf as it s 
detection medium that is expected to be completed in 2011 (llceCube Collaboration! 120071) . 
and KMSNeT, a detector to be built in the Mediter ranean Sea (|Katd!2006l: !Kapp es etaD 
2007! ) with the combined efforts of t he ANTARES (|Katz]!2004l ). NEMolSaEienza 2006,) , 
and NESTOR jAggouras et al.! !2005ll collaborations. Additionally, the AMANDA detector 
( Ahrens et al. 20041 ) . a pathfinder for IceCube covering a volume of 0.015 km^ of Antarc- 
tic ice that was comp leted in 2000, has already produced a number of interesting results 
(|Berghaus et al.l!2007! l. 



2.4.2 Super-K as a neutrino telescope 

The Super-K volume is a large tank of water with a volume of 5 x 10^^ km^. It was 
designed primarily to detect neutrinos from the sun, supernovae, and particle accelerators 
- its volume is quite small compared to the neutrino telescopes currently being built whose 
primary science target is the high-energy astrophysical neutrino flux. Thus the chances of 
actually seeing this signal in Super-K are not very high. 

Nonetheless, it is still an interesting exercise to look for it in the Super-K data: it pro- 
vides an important consistency check, and Super-K can complement IceCube and KMSNeT 
in terms of sky coverage and livetime. Furthermore, observing a high-energy neutrino 
event in Super-K would provide much more detailed directional information than the more 
sparsely instrumented kilometer-scale detectors. 

To this end, a number of neutrino astronomy studies have b een undertaken by the Super- 



K co llaboration, includi ng searches for neu trino point sources (|Abe et al.l!2006! : !Desai et al 



20081 ). supernova bursts ( Ikeda et al. 2007i ). and annihilations of WIMPs in the Sun, Earth, 
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and galactic center ( Desai et al. 20081 ) . Here we focus on one such study: chapter [3] describes 
a search for a diffuse neutrino flux from unresolved sources performed by looking for ex- 
tremely high energy events in 1679.6 live days of Super-K data and comparing the observed 
flux to the expected atmospheric neutrino background. We use Super-Kamiokande's highest 
energy data sample to set an upper limit on the diffuse flux of neutrinos from astrophysical 
sources. 
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Chapter 3 

Search for diffuse astrophysical 
neutrino flux using ultra-high 
energy upward-going muons in 
Super-Kamiokande I 



This chapter is adapted from the paper "Search for diffuse astrophysical neu- 
trino flux using ultra-high energy upward-going muons in Super-Kamiokande 
I" by Molly E. C. Swanson et al. (the Super-Kamiokande Collaboration), 
which was previously published in the Astrophysical Journal 652, pp. 206-215 



(ISwanson et alJl2006l ^ 



3.1 Introduction 



Neutrinos have great potential as astronomical messengers. The detection of neutrinos from 
astronomical sources has led to many significant discoveries. Two sources of extraterres- 
trial neutrino s in the MeV energy r ange have been detected so far: the S un ( Davis 20031 : 
Koshiballiooi ) and supernova 1987a (jffirata et al.lll987l : Eionta et al.lll987l ). The next step 



for neutrino astronomy is to detect neutrinos in the GeV-PeV energy range, which will open 
a new window on the high energy universe. A wide variety of astrophysical phenomena are 
expected to produce extremely high energ y neutrinos, ranging from active galactic nu - 
clei (AGNs) and gamma-ray bursts (GRBs; Halzen and Hooper 20021 : Gaisser et al. 19951 ) 
to more exot ic sources such as dark matter annihilation or decays of topological defects 
(IStastol 1200^ 1. 

The flux of neutrinos at such high energies is quite small; therefore, large-scale detectors 
are required. One effective technique for observing high-energy neutrinos with an under- 
ground detector is to look for muons produced by or i?^ interacting in the surrounding 
rock. (Throughout this chapter, "muons" will refer to both and /i~.) The muon range in 
rock increases with muon energy, which expands the effective interaction volume for high- 
energy events. Downward-going neutrino-induced muons cannot be distinguished from the 
much larger flux of downward cosmic-ray muons, but since cosmic ray muons cannot travel 
through the entire Earth, upward-going muons are almost always neutrino-induced. Thus, 
upward-going muons provide a suitable high-energy neutrino signal. 
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At muon energies above 1 — 10 TeV, the upward-going muon flux due to neutrinos 
from AGNs is expected to exceed the upward-going mu on flux due to atmospheric neutri- 
nos ( Stecker and Salamon 19961 : Mannheim et al. 2001 ). This cosmic neutrino flux could 
be detected either by searching for point sources of high-energy neutrinos or by detecting a 
diffuse, isotropic flux of neutrinos from unresolved astrophysical sources. A diffuse cosmic 
neutrino flux would be observed as an excess to the expected atmospheric neutrino flux 
at high energies. In this analysis, we focus on searching for a diffuse flux of upward-going 
muons due to neutrinos from astrophysical sources using the highest energy data sample in 
Super-Kamiokande (Super-K). 

This study complements other Super-K searches for as trophysical poin t sources of high 
energy neutrinos that use data over a larger energy range ( Abe et al. 20061 ) . In this chapter 
we describe a search for evidence of a high energy astrophysical neutrino flux in Super- 
K's highest energy upward-going muon sample. In §3.21 we describe the Super-Kamiokande 
detector, and in ^3.3l we give the details of how we selected candidate events from Super-K's 
ultra-high-energy sample. We evaluate our selection process with Monte Carlo in §3.41 and 
calculate the observed upward-going muon flux in ^3.51 Section 13.61 and ^3.71 discuss the 
background due to the atmospheric neutrino flux. Based on the results, we set an upper 
limit in ^3.81 and conclude in ^3.91 Any necessary estimates and approximations have been 
made so that they lead to a conservative result for this upper limit. 



3.2 The Super-Kamiokande Detector 

The Super-K detector, shown in Figure 13-11 is a cylindrical 50 kiloton water Cerenkov 
detector, located in the Kamioka-Mozumi mine in Japan. It is 41.4 m tall and 39.3 m 
in diameter. The detector was constructed under the peak of Mount Ikenoyama, which 
provides an average rock overburden of 1000 m (2700 m water equivalent). Its geodetic 
location is at 36.4° north, 137.3° east, and altitude 370 m. 

Super-K consists of two concentric, optically separated detectors. Until 2001 July the 
inner detector (ID) was instrumented with 11,146 inward- facing 50 cm diameter photomul- 
tiplier tubes (PMTs). The outer detector (OD) is a cylindrical shell of water surrounding 
the ID and is instrumented with 1885 outward-facing 20 cm diameter PMTs. Between the 
ID and the OD, there is a 50 cm thick shell. Photons coming from this region will not be 
detected by either the OD or the ID, so we refer to it as the insensitive region. 



More details about the detector can be found in lFukuda et al.l (|2003l ) . The data sample 
used in this analysis was taken from 1996 April to 2001 July, corresponding to 1679.6 days 
of detector livetime. This data run is referred to as SK-I. 

Super-K is primarily designed to detect lower energy neutrinos from the Sun, the at- 
mosphere, and particle accelerators but can potentially detect the extremely high energy 
neutrinos expected from astrophysical sources as well. This chapter focuses on the events 
at the highest energy end of Super-K's detection range. 



3.3 Event Selection 

The ultra-high-energy sample in SK-I consists of events that deposit > 1.75 x 10^ photo- 
electrons (pe) in the ID. In the low-energy regime, on average about 9 pe are recorded by 
the ID PMTs for each MeV of energy deposited in the tank; the electronics for the ID PMTs 
saturate at about 300 pe. Thus an event with > 1.75 x 10^ pe in the ID corresponds to a 
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Figure 3-1 Schematic drawing of the Super-K detector. Inset at bottom right shows the 
location of the detector within Mount Ikenoyama (cutaway view). 
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Figure 3-2 (a) OD-based muon trajectory fit applied to an example MC downward-going 
muon event, (b) OD-based fit applied to an example MC upward-going muon event. The 
size of the circle around each point is proportional to the charge detected in the PMT. 



minimum energy deposition of approximately 200 GeV, but the actual energy deposition 
could be much higher, since the saturation effect prevents all of the produced pe from being 
recorded. 

At high energies, muons have some probability to lose energy through radiative processes 
such as bremsstrahlung, resulting in an electromagnetic shower that deposits large quantities 
of pe in the detector. For comparison, a muon that traverses the maximum path length 
through the ID (50 m) but does not produce any electromagnetic showers will deposit 
approximately 11 GeV via ionization energy loss, corresponding to ~ 10^ pe deposited in 
the ID. Thus a high-pe cutoff offers a means of selecting high energy events. 

At the high-pe threshold of > 1.75 x 10^ pe, the high level of saturation in the ID 
PMT electronics can cause Super-K's precision muon fitting algorithms to fail. Therefore, 
these ex tremely energetic events are not included in other stud ies of upward-going muons 
in SK-I (iDesai et ahlbood : iFukuda et al.]ll999l : lAshie et ahllioosl ). In this study we analyzed 
this ultra-high-energy data sample separately using a different fitting method based on 
information from the OD. 



3.3.1 Outer Detector Linear Fit 

SK-I's ultra-high-energy data sample contains a total of 52214 events. Most of these are 
either very energetic downward-going cosmic-ray muons or multiple muon events where two 
or more downward-going muons hit the detector simultaneously. In order to select candidate 
upward-going muons from this sample, we applied a simple linear fit to the OD data for 
each event. A linear fit was done on the z-position of each OD PMT versus the time it 
fired, weighted by the total charge in the PMT. Example fits of simulated downward-going 
and upward-going muon events are shown in Figure 13-21 The slope of this fitted line is an 
estimate of — cos O, where is the zenith angle of the incoming muon. A positive fitted 
slope indicates that the muon is upward-going. A similar linear fit was done on the x- and 
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y-positions to determine the full muon trajectory through the detector. 

Since this fitting method is based on the OD (which has a lower resolution than the 
ID), it is not as precise as the muon fitting algorithms used in the lower energy upward- 
going muon analysis. However, it works even when the ID PMT electronics are completely 
saturated and the precision ID-based algorithms fail. 

3.3.2 Selection Cuts 

To select candidate upward-going muons, we applied the OD-based fit to all 52214 events 
in the ultra-high-energy data sample. A cut of cosQ < 0.1 was used to eliminate the bulk 
of the downward-going single and multiple muon events. In addition, the fitted trajectory 
was required to have a path length of > 7 m in the ID. 

To ensure high-quality fit results, we looked at the number of OD PMTs hit near the 
projected entry and exit points. For a true throughgoing muon, there will be a cluster of 
hit PMTs around the entry and exit points. If the fit is accurate, the projected path should 
pass through both of these clusters, so we made an additional cut on the number of OD 
PMTs hit within 10 m of the projected OD entry and exit points, A^oDentry and A^oDcxit- 
We required both AoDcntry and AoDcxit to be 10 or greater, indicating significant clusters 
of OD PMTs that were hit near the entry and exit points of the fitted trajectory. 

Events that do not have AoDentry and AoDexit ^ 10 are generally either stopping muons, 
partially contained events, or poorly fitted throughgoing muons. Stopping muons are muons 
that stop in the detector and only form an OD entry cluster. They typically have energies 
of 1-10 GeV, well below the range of the expected astrophysical signal, so it is appropriate 
to discard events that look like stopping muons. Partially-contained events are neutrino 
interactions that take place inside the detector and only form an OD exit cluster — they 
are not part of the upward- going muon flux incident on the detector, so we want to discard 
these as well. This cut also occasionally eliminates inaccurately fitted throughgoing muon 
events, which reduces the efficiency somewhat but improves the accuracy of the fit results. 

Another possible type of event that can masquerade as a throughgoing muon is a 
partially-contained event with multiple exiting particles. Such an event will create two 
(or more) clusters of hit PMTs in the OD, which could be mistaken as the entry and exit 
points of a throughgoing muon. In order to eliminate such events, we looked at the timing 
between the OD and ID entry points. If the event is truly a throughgoing muon, the OD 
PMTs near the entry point should fire before the ID PMTs. We determined OD and ID 
entry clusters via a simple time-based clustering method and evaluated the mean time of 
hits within 6 m of the OD and ID entry points, tioentry and toDentry 

In an ideal measurement, tiDentry < ioDentry would indicate that the perceived entry 
cluster is actually caused by an exiting particle. However, the timing determination is 
complicated by an effect known as prepulsing. Prepulsing occurs when there are so many 
photons incident on the PMT that they are not all converted to photoelectrons at the pho- 
tocathode - some photons hit the first dynode instead and are converted to photoelectrons 
there. When this happens in an ID PMT, the time it was hit appears earlier than it actually 
was, which will artificially reduce iiD entry, making the ID light appear earlier compared to 
the OD light. If this occurs for an ultra-high energy throughgoing muon event, it could 
cause tiDentry — ioDentry to be negative. To allow for this effect, we used a fairly loose cut 
of —40 ns: if tiocntry — ioDcntry < "40 US, indicating very early ID light, then the event was 
rejected cis db likclv neutrino interaction in the ID. 

After these cuts on cos 6, path length, Aooentry, Aooexit, and tioentry - ioDentry were 
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Table 3.1 Visual scan of candidate upward-going muon events 



Visual scan classmcation 


AT 1 J? J. 

JN umber of events 


Multiple muon events 


164 


PMT "flashers" (malfunctioning PMTs) 


5 


Other noise events 


2 


Downgoing muons (manual fit cos > 0) 


171 


Upgoing muons (manual fit cos < 0) 


1 


Total 


343 



applied to the 52214 events in the sample, 343 candidate events remained. These remaining 
events were then evaluated by a visual scan and a manual direction fit by two independent 
researchers to select events with cos < 0. 

The visual scan eliminates events that can pass the automatic reduction but are obvious 
to the trained human eye as noise, i.e., mainly downward-going multiple muon events and 
"flashers." Multiple muon background events occur when two or more downward-going 
muons from an atmospheric shower pass through the detector simultaneously. They have 
extra energy deposition (and thus are expected to be more common in the ultra-high-energy 
sample) and typically give poor CD fit results due to multiple OD clusters, but they are easy 
to identify visually. Flashers are events caused by malfunctioning PMTs that emit light and 
create characteristic patterns. After the multiple muon events and flashers are removed, the 
manual direction fit separates truly upward-going muons from mis-fitted downward-going 
muons. 

From the 343 candidates, only one event passed the visual scan and manual direction 
fit selection as being truly upward-going. The breakdown of the visual scan and manual fit 
classifications is shown in Table 13.11 

This upward-going muon event selected from the > 1.75 x 10^ pe sample is the ultra-high 
energy upward-going muon signal observed by SK-I. This event occurred on 2000 May 12 at 
12:28:07 UT and deposited 1,804,716 pe in the ID. Based on the manual fit results, the path 
length through the ID was 40 m, and the zenith angle was cos = —0.63, corresponding to 
a direction of origin of (R.A.,decl.) = 20''38™, -37°18'. 

3.4 High-Energy Isotropic Monte Carlo 

In order to calculate the observed muon flux, we need to determine the resolution and 
efficiency of the OD-based fit and other cuts on high-energy muons, and we need to estimate 
the probability that a muon of a given energy will deposit > 1.75 x 10*^ pe in the ID. 

To determine these quantities, we generated a high-energy isotropic Monte Carlo (MC) 
sample. This MC consists of an isotropic flux of muons in monoenergetic bins impinging 
on the Super-K detector, representing a flux of muons from neutrino interactions in the 
surrounding rock. The MC simulates the detector response - i.e., the time hit and charge 
deposited for each PMT. 

Seven monoenergetic bins were used, with muon energies ranging from 100 GeV to 
100 TeV, and 10000 events with a path length in the ID of > 7 m were generated in each 
bin. The simulation was performed using a GEANT-based detector simulation. GEANT's muon 
propagation has been shown to agree with theoretical predictions up to muon energies of 
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3.4.1 Resolution and Efficiency of Event Selection 

The muon trajectory fitting algorithm discussed in §3.3.11 was apphed to the high-energy 
isotropic MC, and the fitted values for cos Q, path length, number of OD PMTs hit near the 
entry and exit points, and direction of trajectory were compared with the true MC values. 
These plots are shown in Fig. 13-31 A plot illustrating the angular resolution of the fit is 
shown in Fig. 13-41 The resolution of 5-12° is poor compared to the typical 1° resolution of 
precision ID fitting algorithms, but those algorithms do not work well on these high-energy 
saturated events. 

The efficiency of the upward-going cuts was estimated by considering all of the MC 
events with true values of cos G < and ID path length > 7 m and then determining the 
fraction of these events that pass the selection cuts described in ^3.3.21 This was done using 
isotropic MC events with > 10^ pe since there were relatively few events in the MC with 
> 1.75 X 10^ pe. (Above 10^ pe the efficiency and resolution do not depend strongly on the 
number of ID pe deposited.) Also, as discussed in §3.4.21 only the energy bins in the range 
3.16 — 100 TeV were considered. The efficiency was calculated as a function of muon energy 
and cos 0. 

Statistical uncertaintieson the effi ciency determination were calculated using the Bayesian 
method discussed by Conway ( 20021 ). Systematic uncertainties due to uncertainties in the 



fitted values of the zenith angle, path length, and number of OD PMTs hit near the en- 
try and exit points were calculated using the resolution histograms for cos 0, path length, 
A^OD entry and A'oDexit shown in Fig. 13-31 First, a la acceptance region containing 68.27% of 
the total area was selected for each resolution histogram. Then the cuts on these parameters 
were varied by 1 a in either direction to determine the effect on the efficiency. 

Another source of systematic uncertainty on the efficiency is prepulsing, which is not 
included in the MC simulation and must be estimated separately. To do this, we com- 
pared the results from the high-energy isotropic MC to a sample of 627 ultra-high-energy 
downward-going muon data events with > 1.75 x 10^ pe in the ID selected by a visual scan. 
For many of these data events, the time difference iioentry — ^ODentry is negative, showing 
evidence of prepulsing not seen in the MC sample. 

Histog rams of the time difference tjY) entry — ^OD entry for the data events and the MC 
events are shown in Fig. 13-51 

Note that for the MC events, the histogram is sharply peaked around a time difference 
slightly above zero, and has very few events with negative time differences. In contrast, the 
histogram for the data events has a wider spread, and has several events with negative time 
differences. The difference between the data and the MC distributions is expected to be 
due mainly to prepulsing. 

To adjust for this, we derive a heuristic correction for prepulsing by drawing a random 
value from the distribution defined by the data histogram and another random value from 
the distribution defined by the MC histogram, and then subtracting the data value from 
the MC value. By repeating this process several thousand times, we build up a distribution 
of the expected difference between the data and MC values of tio entry — ^OD entry Then 
for each event in the MC sample, a random value is drawn from the distribution of the 
expected difference between data and MC, and then adding it to tin entry — ^OD entry for the 
MC event. The result is a rough approximation of what the time difference might be if the 
MC simulation included prepulsing. 
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Figure 3-3 Resolution of OD-based fit on determining track parameters for events from the 
high energy isotropic MC with > 10^ pe in the inner detector. Parameters shown are cos 0, 
path length, and the number of OD PMTs hit within 8 meters of the OD entry and exit 
points. The 1 a acceptance regions contain 68.27% of the total area under each histogram 
in each plot. 
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Figure 3-4 Angular resolution of OD-based fit for events from the high energy isotropic MC 
with > 10^ pe in the ID. 



UHE downgoers 

UHE isotropic MC events before heuristic prepulse correction 
UHE isotropic MC events after heuristic prepulse correction 




Figure 3-5 Distribution of the time difference between OD entry and ID entry clusters 
for UHE downgoing muon events from data, compared to the distribution for the UHE 
events in the high-energy isotropic MC. After applying a heuristic correction to account for 
prepulsing, the MC distribution looks similar to the data distribution. 
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The histogram of the heuristically-corrected ultra-high energy MC events is also shown 
in Fig. 13-51 — the distribution looks similar to the data events, as expected. To estimate 
the systematic uncertainty on the efficiency calculation due to prepulsing, we applied this 
adjustment to the events in the MC sample and found that an additional 1% of the MC 
events would be cut after accounting for prepulsing, which lengthens the negative-side error 
bars on the efficiency estimates by approximately 0.01. 

Finally, one additional correction must be made: the efficiency has been estimated at 
different values for E^, but the flux calculation is for the flux above a threshold energy -B™™. 
Since the efficiency decreases with energy, e(£'^,0) is an overestimate for e (> £"^"^,0). 
This does not lead to a conservative upper limit for the flux, so a correction must be made. 
This requires knowledge of the energy spectrum of the expected signal, so we modeled the 
signal as an isotropic flu x of neutrinos with d^^, {E^)/dE,j oc (a plausible astrophysical 
spectrum lGaissej|l990l ). used the method of §3.71 to estimate the muon flux <I>^ (> £'™™), 
and used this muon flux to extrapolate our efficiency calculations to find e (> £™™, ©) • 

The efficiency as a function of energy was estimated using a linear fit on e(S^,0) 
vs. log-E^ in each angular bin, using the MC results from the energy bins in the range 
3.16 — 100 TeV. This gives an estimator e {E^, 0) for the efficiency: 

e{E^,Q) = A{Q) + B{Q)\o^E^, (3.1) 

where A and B are the fit parameters. 

Then e (> S™'"^, G) is estimated by weighting the fitted efficiency e (-E^, G) by the E^^ ^ 
model muon flux. 



e (> Ef"", Q) = , ^ . (3.2) 

This procedure yields small downward adjustments (0.6—9%) to the calculated efficiency 
in each angular bin. 

Additionally, we must also account for the efficiency of the visual scan and manual fit 
procedure. To do this, we did a visual scan of the 605 events from the high-energy isotropic 
MC with true cosG < 0, true ID path length > 7 m, and > 1.75 x 10^ pe in the ID and 
found that 10 of these events were eliminated in the visual scan. This gives an efficiency of 
roughly 98%, so we adjust our efficiency results by a factor of 0.98. 

The final results for the efficiency of our cuts on upward, throughgoing, > 1.75 x 10^ pe 
muons (including all corrections discussed above) are plotted in Figure | 



3.4.2 Ultra-High-Energy Fraction 

In order to set an upper limit on the flux of upward-going muons from cosmic neutrinos 
one must make an inference about the energies of the upward-going muons in the ultra- 
high-energy sample. Above energies of ~ 1 TeV, muon energy loss in water is dominated 
by radiative processes such as bremsstrahlung, so high-energy muons have some probability 
of depositing large numbers of photoelectrons in the Super-K detector and contributing to 
the ultra-high-energy sample. Since this energy loss is not continuous, it is not possible to 
estimate the muon energy for a single ultra-high-energy upward-going muon event. Rather, 
MC is used to make a statistical statement about the energies of the muons that make up 
the > 1.75 X 10^ pe sample. 
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Figure 3-6 Efficiency of our data reduction on upward, througligoing, > 1.75 x 10 pe muons 
as determined by the high-energy isotropic MC. Error bars include statistical and systematic 
errors. 



The high-energy isotropic MC has been used to determine the fraction k (E^) of muons 
with energy Ef, that will deposit > 1.75 x 10^ pe in the ID, thus contributing to the ultra- 
high-energy sample. Results are shown in Tablel3.2[ Stat istical uncertainties were calculated 
using the Bayesian method discussed by IConwavl (j2002l ). As can be seen in Table [3121 the 
three lowest energy bins — from 100 GeV to 1 TeV — do not make a significant contribution 
to the > 1.75 X 10^ pe sample. Hence, the rest of this analysis was done using only the four 
highest energy bins — from 3.16 to 100 TeV. 



These results for k (E^) are used in §3.5l to calculate the flux <I>p 



'> Since A; (^^ 



increases with energy, it is an underestimate for k (> -E™™) , which leads to a conservative 
upper limit for <I)^ (> ii^™™). 

Also, since this MC only simulates the muon and not the actual neutrino interaction, 



Table 3.2 Fraction of high-pe events k (E^) in high energy isotropic MC 





number of MC events 


k{E^) 




(out of 10000) with 


(with statistical 


(TeV) 


> 1.75 X 10^ pe in ID 


uncertainties) 


0.1, 0.316, 1 





0.0000 

-0.0000 


3.16 


35 


0.0035 ± 0.0006 


10 


113 


0113 +°-°°^^ 
^■^^^'^ _o.ooio 


31.6 


293 


0.0293 ± 0.0017 


100 


879 


0879 
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Table 3.3 The flux of ultra-high energy upward-going muons as observed by SK-I. 





(> i^r') 


(TeV) 


(cm ^s "^sr ^) 



3.16 


2.64 X 10" 


-14 


+16.1% 
-17.9% 


10 


8.23 X 10" 


-15 


+9.48% 
-9.73% 


31.6 


3.25 X 10" 


-15 


+6.67% 
-6.40% 


100 


1.10 X 10" 


-15 


+4.96% 
-4.09% 



Note. — These fluxes include both atmospheric background and potential astrophysical signal at each 
threshold energy. 



the effect of lower energy debris from deep inelastic scattering events that make it into 
the detector has been neglected, again leading to an underestimate of k{Efj_). This effect 
is expected to be small in the energy range considered here — above E"^ = 1 TeV at 
the detector entry, over 80% of the neutrino-induced upward-going muons come from over 
200 m away from the detector, so most of the debris is absorbed by the surrounding rock, 
as determined using the atmospheric neutrino MC discussed in ^3.61 



3.5 Flux Calculation 

The flux of upward-going muons above a threshold energy -E™™ is given by 

1 " 1 

where n is the total number of upward-going muon events observed and Qj is the zenith 
angle of the jth event. The efficiency e (> E™™,0j) and the ultra-high-energy fraction 
k (> -E™''^) are calculated in §3.41 T is the detector livetime, which is 1679.6 days for 
SK-I. A (Qj) is the effective area of the Super-K detector perpendicular to the direction of 
incidence for tracks with a path length of > 7 m in the ID. The average effective area of 
the detector is ~ 1200 m^. 



Equation (j3.3p has been applied to the detected upward-going muon event discussed in 
g3X2]to calculate (> E^"^) for E^''^ in the range 3.16 - 100 TeV. Results are shown in 
Table [3^ Systematic uncertainties include a 0.1% uncertainty on the live time T, a 0.3% 
uncertainty on the effective area A, the total efficiency uncertainties shown in Figure 13-61 
and the statistical uncertainties on k shown in Table 13.21 This flux includes both the 
potential signal from astrophysical neutrinos and a background of atmospheric neutrinos. 
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3.6 Expected Atmospheric Background from Monte Carlo 



When searching for neutrinos from astrophysical sources, the dominant background is the 
high-energy tail of the atmospheric neutrino spectrum. Atmospheric neutrinos are pro- 
duced by decays of pions and kaons formed when cosmic rays interact with particles in the 
atmosphere. In order to set a limit on the flux of cosmic neutrinos, the expected flux of 
upward-going muons due to the atmospheric background must be carefully analyzed. To do 
this, we have used an atmospheric neutrino MC that is a 1 00 yr equiva lent sample of events 
due to the atmospheric neutrino flux. The neutrino flux in lHonda et al . (2004) was used up 
to neutrino energies of 1 TeV. At 1 TeV, the calculated flux in IVolkoval dlMfll) was rescaled 
to the Honda et al. flux. Above 1 TeV, the rescaled flux from Volkova was used up to 100 
TeV. Neutrino int eractions were modeled using the GRV94 parton distribution functions 
(jGluck et al.lll995l ). and muon propagation through the rock and water was modeled usiii g 
GEANT. Further details on the atmospheric neutrino MC can be found in lAshie et al.l (|2005l ). 
No correction is m ade for neutrino o scillations, because based on the oscillation parame- 
ters determined in lAshie et al.l (|2005l ^. the neutrino oscillation probability is negligible for 
neutrinos above 1 TeV. 

The model flux is weighted to represent 3 years of solar minimum, 1 year of changing 
activity, and 1 year of solar maximum, which corresponds to the actual solar activity during 
the run of SK-I. This atmospheric MC is split into two parts: a partially-contained/fully- 
contained (PC/FC) sample, which consists of events with neutrino interaction points inside 
the ID plus a shell 50 cm thick surrounding the ID (the insensitive region) , and an upward- 
going muon sample, which consists of events with neutrino interaction points outside the 
ID. Note that these two samples overlap because they both cover the 50 cm insensitive 
region. 

The OD-based fit was applied to the events in the atmospheric MC, using the same cuts 
that were applied to the SK-I data. A total of 11 MC events passed the > 1.75 x 10^ pe, 
cos© < 0.1, path length > 7 m, A'^oDentry and A^oDexit ^ 10, OD/ID timing, and manual fit 
cuts. Out of these 11, 2 are from the PC/FC sample, both with interaction points inside 
the ID. The remaining 9 events are from the upward-going muon sample: 3 events with 
interaction points in the 50 cm insensitive region, 1 event in the water of the OD, and 5 
events in the rock surrounding the detector. 

All of these background events are deep inelastic scattering (DIS) events where an 
interaction between a muon neutrino and a nucleon produces a muon plus a spray of lower 
energy particles. The 6 events with interaction points within the detector (ID or OD) have 
muon energies of 0.1 — 0.8 TeV, and the 5 events occurring in the rock have muon energies 
of 1 — 20 TeV. This difference in the energy range can be understood as follows: For DIS 
events occurring a long distance (> 2 m or so) from the detector, only the muon will reach 
the detector since the lower energy debris will be absorbed by the surrounding rock, but 
for nearby events or events occurring in the water of the OD, some of these lower energy 
particles will enter the detector as well. This means that nearby events can be included in 
the > 1.75 X 10^ pe sample with lower muon energies than more distant events. 

Since the insensitive region is covered by both the PC/FC and the upward-going muon 
MC samples, we divided the 3 events originating from this region in half, for a total of 1.5 
events in the insensitive region. This gives a total of 9.5 MC events in 100 yr of simulated 
live time. Scaling the 100 yr MC to SK-I's live time of 1679.6 days gives an expected 
background of 0.44 events due to atmospheric neutrinos during the operation of SK-I. 
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Table 3.4 Systematic uncertainties in atmospheric neutrino background 



Source of Uncertainty 


Uncertainty 


Statistical 


31% 


Absolute normalization of atm v flux 


10% 


Primary spectral index 


37% 


Neutrino cross-section uncertainty 


10% 


Total uncertainty in background flux 


50% 



The statistical uncertainty in this background measurement of the MC events is 

= 31%. (3.4) 



1 

9^ 




There are also significant systematic uncertainties: the normalization of the atmospheric 
neutrino flux has a theoretical uncertainty of ±10% at neutrino energies below 10 GeV 
(lAshie et al.ll2005l ). M order to extend this to the energy range of the expected background, 



we must also account for the uncertainty of 0.05 in the spectral index of the primary cosmic 
ray spectrum above 100 GeV, which leads to a . 05 un certainty in the spectral index for 



atmospheric neutrinos above 10 GeV (jAshie et al.ll2005l ). 



To determine how much this uncertainty affects our result for the background, we con- 
sider the spectral index 7 to be known at EY"°^ = 10 GeV, and we calculate the uncertainty 
of the total flux above a threshold energy i?™™ = 10.6 TeV, the average neutrino energy 
of the 100 year MC events passing our cuts. For a differential flux of d(^y/dEy oc E,y"' , the 
uncertainty in due to the uncertainty in 7 is given by 



With 7 = 3.7, (57 = 0.05, = 10.6 TeV, and ^J^™* = 10 GeV, this formula yields 

b^yj^v = 0.37. Thus the primary spectral index uncertainty gives us an additional ±37% 
uncertainty on the atmospheric neutrino flux 

Finally, the neutrino cross-section at high energies is thought to be known to within 10% 
or less, so we include an additional 10% uncertainty to account for this. These uncertainties 
are summarized in Table [33] and lead to a total uncertainty on the background of 50%. 

The atmospheric MC does not include neutrino oscillations, so this is another potential 
systematic uncertainty. However, for neutrinos above 1 TeV, the neutrino oscillation prob- 
ability is negligible, so it does not make a significant contribution for this analysis. Other 
potential errors (uncertainties in the simulation of the SK detector and varying hadron 
multiplicities in different deep inelastic scattering models) were tested and shown to not 
make a significant contribution to the systematic uncertainty. 

Another potential background source of high-energy neutrinos not included in the 100 
yr atmospheric MC is the prompt atmospheric neutrino flux, which arises from decays 
of short-lived charmed particles produced when cosmic rays interact with particles in the 
atmosphere. This flux is not as well-understood as the conventional atmospheric neutrino 
flux from decays of pions and kaons, but it is expected to have a harder spectrum and 
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therefore is expected to become more important as we push towards higher energy scales. 

Based on the 100 yr atmospheric MC, we calculate an expected background for this 
analysis of 0.44 it 0.22 events, compared to the 1 event observed. However, there are three 
effects that this MC does not take into account: it does not include neutrinos over 100 TeV, 
it does not include the prompt atmospheric neutrino flux, and it does not account for 
attenuation of neutrinos passing through the Earth. We account for these issues by making 
corrections based on an analytical calculation discussed in ^3.71 

Also, it is important to note that the atmospheric background — both conventional and 
prompt — comes from a lower energy range than that for which we expect to observe the 
possible signal of neutrinos from astrophysical sources. Roughly speaking, the peak 90% 
of the expected upward-going muon events come from muons with energies in the range 
Ef^ = 0.02 — 10 TeV for conventional atmospheric neutrinos, and = 0.2 — 200 TeV for 
prompt neutrinos. As we learned from the 11 background events in the 100 yr atmospheric 
MC, muons with energies E^ < 1 TeV contribute to the 1.75 x 10^ pe sample mainly via 
debris from DIS events very close to the detector rather than from catastrophic energy loss 
of the muon. 

Since there are many more low-energy events in the atmospheric spectrum, they will 
dominate even though each one only has a tiny probability of depositing a large amount 
of energy in the detector. In contrast, the peak 90% of the expected muon events from a 
harder spectrum — a hypothetical E~'^ astrophysical flux — come from the range E"^ = 
7 — 6000 TeV. Thus, even though the atmospheric flux is very small in the energy range 
Efj^ = 3 — 100 TeV and above where we are setting our limit, the high-pe tails of the 
distribution from lower energy events dominate our background simply because of the much 
larger flux of lower energy atmospheric events. 



3.7 Analytical Estimate of Expected Muon Flux 
3.7.1 Method for Calculating Muon Flux 

In order to better understand our observed flux of high-energy upward-going muons, we 
developed a method to calculate the expected upward-going muon event rate due to a 
predicted flux of neutrinos. We have used this to plot curves for theoretical muon fluxes in 
Figure 13-91 and also to adjust the atmospheric background calculated with MC in §3.61 by 
correcting for effects not included in the simulation. 

To convert a model neutrino flux into an expecte d up ward-going mu o n eve nt rate, we 
follow the calculation detailed in ICaisser et all \l99^ ) and ICandhi et al.l (|l996l l. The flux 



of muons (> i^"™™) above an energy threshold -E™™ is given by 

% (> i?r) = / dE^P^ {E,, Ef-) ^tM, (3.6) 

where {Ey, -E™™) is the probability that an incoming neutrino with energy Ey will 
produce a muon with energy above the threshold E'™" at the detector, and d<^^ {Eu)/dEy 
is the differential neutrino flux averaged over solid angle and reduced by an exponential 
factor due to attenuation of the neutrinos as they pass through the Earth. 
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given by 



dndE^ 



(3.7) 



Here, the integral over solid angle extends over cos < 0, [d^<I>jy {Ey, 0, (/>)] / [dQdEy] is the 
angle-dependent neutrino flux, A (0) is the effective area of the SK detector at a zenith angle 
, and is the average effective area j^-^ dcos (0) A (0). The exponential factor accounts 
for attenuation of neutrinos as they pass through the Earth, which becomes important at 
higher energies. 

Na = 6.022 X 10^^ mol~^ =6.022 x 10^^ cm~'^ (water equivalent) is Avogadro's number, 
a [Ey) is the total neutrino-nucleon cross-section, and z (0) is the colu nin depth through the 
Earth at a zenith angle 0. The col umn depth has be en calculated in iGandhi et al. 
using the Preliminary Earth Model ( Dziewonski 19891 ) — we adopt this column depth here. 



As discussed in Gandhi et al. 



there is an uncertainty in choosing a value for 
the cross-section a {Ey) to use in the neutrino attenuation factor. Neutrinos interacting in 
the Earth via charged-current (CC) interactions will disappear from the incoming neutrino 
beam. However, neutrinos interacting via neutral current (NC) interactions will only be 
scattered and reduced in energy, so they may still contribute to the overall neutrino flux. 
Thus, we can put upper and lower limits on the appropriate value of a to use in Eq. 13.71 
ccc < c < CTcc + o"7vc. Here we calculate the angle- averaged flux using both the upper 
and lower limits to estimate the resulting uncertainty, and then average the resulting fluxes 
together. 

If the model neutrino flux is isotropic, it can be taken out of the angular integral in 
equation ()3.7p . and we can derive a shadow factor S (E^) given by 



(3. 



S (E^) = f dcos (0) A (0) e(-^(e)-(^^)^A). 
Thus, for an isotropic neutrino flux, equation (j3.6p simplifies to 



(> = / 

JE] 



dEyP^ iyEy, E"™" 



)S{E, 



dE^ 



(3. 



The probability {E,y, -E"™™) depends on both the charged-current neutrino cross- 
section and the energy lost by the resulting muon as it propagates through the Earth. It is 
given by 

dace 



[Eu, E^ 



Na 



E, 



dE^j. {E^, Eu) R {E^, E™") 



(3.10) 



where dace/ dE^ is the differential (with respect to muon energy) charged-current neutrino 
cross section and R {Efj_, i?™'") is the average range that a muon with energy E^ will travel 
before its energy is reduced to -E™™. 

In our calculation, we ca lculate the cross-section using t he GRV94 parton di stribution 
functions ( Gluck et al. 1995 ) applied with the code used in Gandhi et al. ( 19961 ) provided 
by M. Reno (2005, private communication), and we use an effective muon range calculated 
using MC methods ( Lipari and StanevI 1991 ). 

We integrated equation (j3.6p numerically using various predicted models of the neutrino 



flux d^^y/dEy to calculate theoretical <I>^ (> E™™) curves. These are shown in Figure 
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for comparison to our experimental limits. 



3.7.2 Analytical Estimate of Expected Background 

Of particular interest here are the model fluxes of the background due to atmospheric 
neutrinos, both conventional and prompt. We use the method discussed in ^3.7.1l to correct 
for the omissions in the 100 yr atmospheric MC discussed in §3.6l by accounting for neutrinos 
over 100 TeV, attenuation of neutrinos in the Earth, and the flux of prompt neutrinos. 

To estimate the number of expected background events with our analytical calculation, 
we need to factor in the number of muons at each energy that will make it above the 
> 1.75 X 10^ pe cutoff. The expected number of events seen by Super-K in livetime T is 
given by 

N = 27rTA,^J^ dEf^ dE^in'^ ^ ^ (3-11) 

where d<I>^ (> E™™) /dE^^^^ is the derivative of the curve calculated by the method in §3.7.11 
and k (.E™"^) is the fraction of muons with energy above E™™ that will deposit > 1.75 x 
10 pe in the ID. We estimated k (E™''') using the results of the high-energy isotropic MC 
discussed in ^3.4.21 The four non-zero values of k (E^) shown in Table 13.21 are used to 

find an estimator /cpL by fitting with a simple power law: log /cpl {E^) = A + B\og [E^] 
. Additionally, we impose a physically motivated high-energy cutoff: since k represents a 
probability, it should never exceed 1. Thus our estimator for k (E™*"') is given by 

k (E--) = max (fepL (E™") , l) (3.12) 

The energy at which kpi^ crosses 1 is approximately 1.5 x 10^ GeV, well above the energy 
range that produces most of the flux. Figure [3^ compares this power-law fit to results from 
100 yr atmospheric MC, illustrating that the power law is reasonable in the energy range 
we are considering. 

We analytically estimated the expected number of events in the 100 yr atmospheric 
MC by starting with the same input neutrino flux, applying equation (j3.6p without the 
exponential neutrino attenuation factor and integrating up to a maximum neutrino energy 
of EJJ^'^^ = 100 TeV. 

As a cross-check on the analytical calculation, we compared the resulting angle-averaged 
muon flux with the data from the 100 year atmospheric MC. Figure 13-81 shows that the 
analytical calculation predicts the shape of the flux curve accurately, but has a slightly 
lower normalization. Since we only use the analytical calculation to derive relative scalings, 
we consider it to be in good agreement with the MC. 

We then used this muon flux in equation (j3.1ip to find the expected number of back- 
ground events in 100 yr and obtained A''ioo = 7.39 expected events. (The 100 subscript 
denotes 100 yr of exposure.) This matches our results from the 100 yr atmospheric MC to 
within statistical uncertainties. 

To estimate the effects of neutrinos above 100 TeV, we repeated the above calcuation 
without the E™'^^ = 100 TeV cutoff and obtained A''ioo = 7.47 events, an increase of 1.1%. 
Including the neutrino attenuation factor as well, we obtained = 7.32 events for a = acc 
and = 7.27 events foru = acc + (^NC-, corresponding to an average of A^ioo = 7.29 ± 
0.02 events, a decrease of 2.4% ± 0.3%. 
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Figure 3-7 Fraction of events k (E^) that deposit > 1.75 x 10® pe in the ID for the 100 
yr atmospheric MC {open circles), the high-energy isotropic MC {filled squares), and a 
power-law fit to the results from high-energy isotropic MC. The power law is a reasonable 
approximation for the dominant energy range of the atmospheric flux {Ef^ = 100 GeV — 
30 TeV). 
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Figure 3-8 The number of upward-going muon events above energy threshold induced 
by atmospheric neutrinos. Predictions from our analytical calculation are compared with 
the results from the 100 year atmospheric MC. 
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Table 3.5 Confidence intervals for the upward-going muon flux due to neutrinos from astro- 
physical sources. 



(TeV) 


9 


0% C.L. ran 


ge 


( 


cm^^s^^sr^ 


') 


3.16 





- 1.03 X 10" 


-13 


10 





- 3.19 X 10- 


-14 


31.6 





- 1.26 X 10- 


-14 


100 





- 4.28 X 10- 


-15 



We also used this calculation to correct for the prompt atmospheric neutrino flux. To 
account for the theoretical uncertainty in the prompt flux due to differences between various 
flux models, we defined a high m odel and a low mode l for the prompt flux that bracket 
the models shown in Figure 1 of iGelmini et al.1 jiool) that are not ruled out by 
imental limits. The models used are discussed in more detail in Thunman et al. 



exper - 
il996l ). 



Zas et al.1 (fl993|) . iRvazhskava et al.1 (120021 ) . iBugaev et al.l |l989l ;) . IPasguah et al.1 (|l999l ). and 



Gelmini et al.l (2000a 



3). Our analytical calculation gives A^ioo = 0.033 events for the low 



model and A'^ioo = 0.94 events for the high model, corresponding to a prompt flux that is 
6.7% lb 6.2% of the flux of conventional atmospheric neutrinos. 

Based on these results, we correct the background estimate made using the 100 yr 
atmospheric MC in ^3.6l bv applying these relative scalings to the result of 0.44 events in the 
SK-I livetime from §3.61 This gives us a final result for the background of 0.46 it 0.23 events 
for the SK-I exposure. 



3.8 Upper Limit for Muon Flux from Cosmic Neutrinos 



Using the observed ultra-high-energy upward-going muon signal of 1 event and the expected 
atmospheric neutrino background of 0.46 it 0.23 events, we have calculated 90% confidence 
upper limits for the upward-going muon flux in the 3.16 — 100 TeV range due to neutrinos 
from astrophysical sources (or any othe r non-atmospheric sources). 

This was done using the method of iFeldman and CousinsI (119981). w i th the systematic 
uncertainti es incorporated usin g the method of Cousins and Highland! ( 1992 ). as imple- 



HiU 



200 



This method incorporates 



mented by Conrad et al. ( 20031 ) and improved by 
both uncertainties in the background flux and uncertainties in the flux factor / relating 
the observed number of events n to the observed flux <I>: ^ = fn. The uncertainty in / 
includes systematic errors in the livetime, effective area, efficiency, and ultra-high-energy 
fraction. For the confidence interval calculation, the largest percent error for each energy 
bin from Table [331 was used as the percent uncertainty in /. The uncertainties in both the 
background and the fiux factor were assumed to have a Gaussian distribution. 

The final results are shown in Tabl e 13.51 and plotted in Figure 13-91 along with models 
of various p ossible signals from AGNs (Stecker and Salamon 19961 : Mannheim et al. 200ll ) 



and GRBs ( Waxman and Bahcalll 1999 ), a s well as the backgroun ds due to atniosphe ric 
neutrinos ( Honda et al. 20041 : IVolkoval 1980l ) and prompt neutrinos ( Gelmini et al. 20031 ). 

The upper limits calculated here are consistent with the models of astrophysical sig- 
nals. Also shown are l imits on a hypotheti cal E~'^ isotropic neutr ino fiux set by MACRO 
(jAmbrosio et al.l boosi ) and AMANDA-II (jAchterberg et al.l lioOTl ). The model neutrino 
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Figure 3-9 Upper limits from this analysis on muon (//^ + ) flux above energy threshold 
S™", compared to various model fl uxes. Models shown for mu on flux du e to astrophysical 



neutr inos are AGN models from SS (IStecker and Salamonlll996l) a nd MPR (jMannheim et al 



l2nnih . and a GRB model from WB (IWaxman and Bahcalllll999l ) Also shown is t he atmo 
spheri c background, as modeled by Honda et alj ( 2004 ) below Ey = 1 TeV and by Volkova 
(|l980l ) rescaled to match the Honda model above Ey = 1 TeV. The upper edge of the 
atmospheric band represents the horizontal flux, and the lower edge represents the vertical 
flux. The background due to muons from prompt atmospheric neu trinos (assumed to b e 
isotropic) is shown for a range of possible models as summarized in iGelmini et aD \2mi \. 
Finally, we also show limits on a hypothe tical E~'^ isotropic neu trino flux set by MACRO 
(jAmbrosio et al.ll2003l ) and AMANDA-H (jAchterberg et al.ll2007l ). The models of the neu- 
trino flux have b een converted to a muon flux with equation (j3.6p usin g the GRV94 par- 
ton di stributions ( Gluck et al. 19951 ) and the effective muon range from Lipari and Stanev 
(|l99ll ). 
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Table 3.6 Approximate upper limits from SK-I on astrophysical neutrinos (f^ + z^^) 



771mm 

% 


90% C.L. upper limit 


Neutrino energy range 


(TeV) 


(GeV cm~^s~^sr~^) 


(GeV) 


3.16 


6.0 X 10"^ 


6.3 X 10^ - 1.4 X 10*^ 


10 


3.7 X 10"^ 


1.7 X lO'' - 2.4 X 10^ 


31.6 


3.2 X 10-5 


4.7 X 10^ - 5.1 X 10^ 


100 


2.6 X 10-5 


1.4 X 10^ - 1.1 X lO'^ 



Note. — Upper limits on [d^^ /dE^). Note that converting from a muon flux limit to a neutrino flux 
limit requires additional assumptions — our limits on the muon flux are shown in Table l375l 



fluxes were converted into muon fluxes using equation (|3.6p as discussed in ^ 

To facilitate easier comparison with other experiments, we also convert our limits on 
the muon flux into approximate limits on the neutrino flux. In order to do this, we assume 
a model neutrino flux that is isotropic and proportional to E~'^. To get an approximate 
neutrino limit, we find normalization factors for an i?"^ muon fiux curve such that the curve 
passes through each of our four limit points in Figure 13-91 and we use these factors to find 
the implied limits on E'j'^ fiux. 

In order to determine the approximate neutrino energy range in which these limits are 
valid, we use equation (13. 6p to determine the neutrino energy range that produces the bulk 
of the muon signal for a given value of the muon energy threshold -E™™- We define the 
energy range as the range that (1) produces 90% of the muon fiux above -E™™ and (2) has 
a higher value of the integrand of equation (j3.6p within the range than anywhere outside 
the range. T his defini t ion is based on the definition of highest posterior density intervals as 
described by Conwav ( 20021 ). 



Explicitly, the energy range is {Ea, Eb), where 



fp " dE^I (E^) 

.9 = ^ \ ' (3.13) 

j£;mill dEyl {Ey) 



with the integrand defined as 



I [E,] = P (Ey, Ef'') S (Ey) (3.14) 



dEy 



as in equation ()3.9p . with 



d£„ oci.-. (3.15) 

Since this requirement alone does not define a unique range, we also require that the value 
of / (E'l/) at any E^ G {Ea, Eh) be greater than the value of I (E^) at any Ei, ^ (£'„, E^). 

The results of these approximations for neutrino limits and energy ranges are shown in 
Table \3M and are also plotted in Figure [3-10^ with the same models and experimental limits 
as shown in Figure l3-9[ To draw our limits on this plot, we chose to draw only the most 
sensitive limit in each energy range. 
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Figure 3-10 Approximate upper limits from SK-I on astrophysical neutrinos (f^ + z^^). Mod- 
els shown are the same as in Figure [3^ Note that converting from a muon flux limit to a 
neutrino flux limit requires additional assumptions - our limits on the muon flux are shown 
in Figure l3-9[ 
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Note that these results for the hmits on the neutrino flux are only approximations made 
to facilitate comparison to other experiments — our primary results are the limits on the 
muon flux shown in Table 13.51 and Figure l3-9[ 



3.9 Conclusions 



In conclusion, we have developed a method for analyzing Super-K's highest energy data to 
search for evidence of high-energy neutrino flux from astrophysical sources. We have done 
a thorough study of the efficiency and the expected backgrounds from this method and 
applied our method to the SK-I data sample. Our study of the highest energy events in 
SK-I does not show evidence of a high-energy cosmic neutrino signal. 

We have set upper limits on the muon flux due t o cosmic neutr ino sou rces. These limits 



are c onsistent with the results of other experiments (jAmbrosio et al . 2003: .Achterberg et al 



20071 ) . It is possible that an astrophysical neutrino signa ,l could be within the grasp o f the 
next gener ation of neutrino detectors suc h as IceCube ( IceCube Collaboration 200?! ) and 
KM3NeT dKatzl I2OO6I : IXaPDes et al.ll2007l ). 
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Part II 

Particle Physics in the Sky 
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Chapter 4 

Doing cosmology with galaxies 



4.1 Cosmology Basics 
4.1.1 What is cosmology? 

Cosmology is the study of our universe as a whole, and thus focuses on the largest scales 
accessible to science. It strives to answer a number of "big questions" : What is our universe 
made of? How did it begin? How did the matter assemble into the structures we see today? 
What is its ultimate fate? In order to address these questions, cosmologists use a wide 
variety of astronomical observations and draw on theory from across essentially all fields of 
physics, ranging from general relativity to quantum field theory. 

One remarkable feature of the study of cosmology is that it provides insights into particle 
physics - i.e., the study of the smallest scales in physics - that could never be observed in 
a terrestrial laboratory. The tremendous energy in the hot, early universe shortly after the 
Big Bang far exceeds the energies that could ever be produced by the most sophisticated 
particle accelerators on Earth. Now that cosmology is becoming a precision science, particle 
physicists are turning to cosmology as a tool, using the study of the largest scales in physics 
to learn more about the smallest. 



4.1.2 The standard cosmological model 

Over the past few decades, c osmology has ev olved from a highly speculative field dubbed "A 
search for two numbers" by Sandage ( 197d ) into a full-fledged observational science, with 
vast quantities of data supporting a detailed theoretical model that is now well-established 
enough to be known as the "standard cosmological model". We are commonly said to be 
living in the age of precision cosmology. However, many aspects of the standard cosmological 
model are quite surprising, and have led to more questions than answers. 
The basic premises of the standard cosmological model are as follows: 



• Our universe is expanding - the gravitationally bound structures in our universe (e.g. 
clusters of galaxies) are all moving away from each other. 

• Our universe used to be much hotter, denser, and smoother than it is today - the 
early universe was a hot soup of quarks and elementary particles. 

• The large-scale structures in our universe grew through gravitational instabilities 
seeded by quantum fluctuations in the early universe. 
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• These quantum fluctuations were stretched to macroscopic size during a phase of 
rapidly accelerating expansion of the early universe called inflation. 

• Ordinary matter (protons and neutrons, a.k.a. baryons) only make up 4% of the mass 
density in our universe. 

• About 21% is made of dark matter, a mysterious substance that is gravitationally 
attractive but does not interact with light. 

• About 75% is made of dark energy, an even more mysterious substance that eff'ectively 
gives rise to a repulsive gravitational force and is causing the current expansion of our 
universe to accelerate. 

This basic picture is consistent with a wide variety of different types of measurements, 
ranging from the tiny fluctuations in the microwave radiation produced by the early universe 
to the rates at which distant supernovae are moving away from us to the clustering patterns 
observed in the distribution of galaxies today. The agreement between all of these different 
measurements is quite remarkable, and has forced scientists to take seriously these seemingly 
preposterous concepts of dark matter, dark energy, and inflation. 

Today the study of cosmology is focused on gaining a deeper understanding of the 
physics behind these concepts. What triggered inflation and how did it stop? What is the 
nature of dark matter? How much of the dark matter might be neutrinos? What sort of 
substance could have such bizarre properties as dark energy? Is dark energy a substance 
at all, or is the observed acceleration really due to a breakdown of general relativity at 
cosmological distance scales? These are the types of questions currently being investigated. 



4.1.3 Zeroth order cosmology 

As a first approximation, our universe can be treated as being homogenous and isotropic - 
that is, the matter in our universe is distributed uniformly and it looks the same in every 
direction. Looking around our neighborhood might lead us to believe that this is not a very 
good assumption: matter is packed together much more densely within the Earth than a 
thousand kilometers above its surface, for example. 

However, cosmologists think on much larger scales: if we average over scales of several 
hundred Mpc (i.e. about a billion light years - that's 100 trillion times the distance from 
the Earth to the Sun), we will smooth out structures as small as a solar system, a galaxy, or 
even a galaxy cluster, and the matter in our universe begins to look quite uniform. Histori- 
cally, the homogeneity and isotropy of our universe on large scales has been a fundamental 
assumption of cosmology, but in recent y ears our surveys of our universe have gotten large 



enough to confirm this directly (see, e.g., Yadav et al. 20051 ). 



The physics of a homogenous and isotropic universe can be described in Einstein's theory 
of general relativity (or more generally, any metric theory of gravity) using the Friedmann- 
Robertson- Walker (FRW) spacetime metric: 



-c^dt^ + a (t) 



dr 



l-Kr"^ 



+ {de"^ + sin^ ed(t)'^) 



(4.1) 



This resembles the more familiar Minkowski metric used in special relativity, given by 



ds' 



-c^dt^ + dr'^ + 



{de^ + sin^ edcf)'^) 



(4.2) 
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Universe with positive Universe with negative cun/ature. 

curvature. Diverging line Lines diverge at ever increasing angles, 

converge at great distances. Triangle angles add to less than 1 80°. 
Triangle angles add to more 
than 180°. 



Universe with no curvature. Lines diverge at 
constant angle. Triangle angles add to 180°. 

Figure 4-1 The curvature of space illustrated for two-dimensional universes. Figure copy- 
right Nick Strobel, p.ttp : //www . astronomynotes . com. Used with permission. 



with two additional pieces: a constant factor K describing the global curvature of space, 
and an arbitrary function of time a (t). The curvature of space is defined by the behavior 
of parallel lines: in flat space {K = 0), parallel lines always remain parallel, in a space with 
negative curvature {K < ), they diverge, and in a space with positive curvature {K > 0), 
they converge. An illustration of this is shown for two-dimensional spaces in Fig. 14- our 
universe is a three-dimensional analog of this. The function a (t) is known as the scale 
factor, which allows space itself to expand or contract as our universe evolves. Here we 
normalize a (t) such that a (to) = 1 at the present time to. Doing cosmology at zeroth order 
fundamentally boils down to determining K and a (t). 

By applying the E instein field equ ations for a homogeneous, isotropic universe to the 
FRW metric (see, e.g., IPodelsonllicm ;). K and a (t) can be linked to to the composition of 



our universe via the Friedmann equation: 



H^^{-] =-_p(t)-_^. (4.3) 



We will discuss each term in this equation in turn. The factor a/a is the expansion 
rate of our universe. It is called the Hubble parameter and denoted by H. The value 
of H today is known as the Hubble constant Hq = a (to) /a (to). This is the parameter 
that Edwin Hubble measure d in his initial observations of the expansion of our universe 
( Hubble and Humason 193ll ). He observed that the further away a galaxy was, the faster 
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it appeared to be moving away from us: 



V = Hod (4.4) 

where f is a galaxy's recession velocity and d is its distance from us. Bubble's initial mea- 
surements gave Hq ~ 550 km s~^Mpc~^; the modern accepted value is Hq = 73 km s~^Mpc~^. 
It is often parameterized using the Hubble parameter h, where h = Hq/ [lOOkm s~^Mpc~^] . 
Cosmological distances are conventionally written in units of /i^^Mpc - we use this conven- 
tion in chapter [6l 

The recession velocity v is inferred from the wavelength shift of the galaxy's spectral 
lines due to the Doppler effect: 

V Aobs — Arest / a r\ 

- ^ = z. (4.5) 

where c is the speed of light, Arest is the wavelength of the spectral line in the rest frame 
of the emitter, Aobs is the wavelength observed at Earth, and the quantity z is called the 
"redshift" since objects moving away from us appear redder. The approximation holds 
when u <C c and the galaxy is close enough that it is sensible to interpret the redshift as a 
recession velocity. However, it is more accurate to think of the redshift as being caused by 
the expansion of space as the light is propagating through our universe. The wavelength 
of the light will get stretched along with space, and thus the wavelength change can be 
directly related to the scale factor a: 

Arest a (to) 

where te is the time the light was emitted from the galaxy. Light from galaxies that are 
further away will take longer to reach us, so it will have a longer amount of time in which 
to get stretched and thus appears more redshifted. Using equations (j4.5p and (j4.6p along 
with the normalization convention a (to) = 1, we can relate z to a: 



l + z 



(4.7) 



Moving on to the first term on the left-hand side, we have the the matter-energy density 
of our universe p{t). The time evolution of this quantity depends on the type of matter. 
For the case of a perfect isotropic fluid with density p and pressure P, the Einstein field 
equations and the FRW metric tell us that (see, e.g., Dodelson 20031 ) 



p + 3-(p+^]=0, (4.8) 

Determining the evolution of p requires a relation between P and p, known as the equation 
of state. Here we will consider three different types of contributions to p: 

• Cold, pressureless matter p^ (includes baryons and dark matter): P = 0, so equa- 
tion (j4.8p yields pm oc . Physically this means that the number of particles stays 
constant while the volume expands as length^. 

• Radiation or relativistic matter p^ (includes photons and neutrinos): P = p/3, so 
equation (|4.8p yields pr oc a~^. Again the number of particles stays constant and 



60 



the volume expands as length^, but we gain an additional factor of since the 
energy of each particle is ^ oc A^^, so the energy per particle decreases as as the 
wavelength gets redshifted. 

• Unknown equation of state, i.e. dark energy P = wp, so equation (j4.8p yields 
PA oc a~^^^~^'^\ Current observations of dark energy are consistent with a cosmological 
constant, i.e. an energy density that remains constant as our universe expands. This 
corresponds to pA oc a*^, so u; = — 1. In the above examples, pressureless matter has 
w = and radiation has w = 1/3. More general parameterizations of the dark energy 
equation of state also allow w to vary in time. 

Finally, the second term on the right-hand side of equation (j4.3p . which includes the pa- 
rameter K defining the overall curvature of space in our universe. One important insight to 
be gained from this equation is that the global geometry of our universe is directly linked 
to its matter content. There is a special value of the matter density that corresponds to 
flat space (i.e. K = 0) that is known as the "critical density" pcr- By plugging K = into 
equation (j4.3p we find that the present value of the critical density is given by 

.c,„ = 5§. (4.9) 

The Friedmann equation takes on a particularly elegant form if we express all densities as 
fractions of the critical density. We define 

n^^^,nr^fi^,nA^^, (4.10) 

PcrO PcrO PcrO 

where again the subscript denotes the present-day value. Now the Friedmann equation 
evaluated at the present is 

Kc^ 

1 = 0,jYi + fij, + Oa ttT' 

Hq 

Using this we can define an effective Q for the curvature: 

Kc 



VlK^l-^m-^r-^k = TpT- (^-H) 



Thus the curvature is determined by the total matter density: if $7^ + ^r + = 1) i-e. the 
total density of matter is exactly equal to the critical density, then Q.^ = and our universe 
is flat. If the total density is less than the critical density, then K <Q and our universe has 
negative curvature. If the total density is greater than the critical density, then K > Q and 
our universe has positive. 

To write our elegant Friedmann equation, we use the definitions in equations (j4.10p and 
(14. lip and the scalings with a of different types of matter density determined by equa- 
tion ()4.8p to re-express equation (j4.3p as 

H =Hl i^^ra-^ + O^a"^ + ^Ka-'^ + J^Aa^^^^^"'^) (4.12) 

This relatively simple differential equation for the scale factor a (t) is the cornerstone of 
zeroth order cosmology. Given values for the parameters Hq, ^r, ^A, and w, one can 
completely describe the evolution of a homogenous universe containing matter, radiation. 
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and dark energy. (Or to be more precise, a dark-energy-like component with an equation 
of state that can be parameterized by a constant w.) 



4.1.4 First order cosmology 

The next step in cosmology is to study the evolution of small deviations from homogeneity. 
The early universe was extremely smooth, with tiny perturbations in the density on the order 
of one part in 10^. As our universe evolved, these small fluctuations grew under the influence 
of gravity into the complex cosmic web we see today. The behavior of these fluctuations is 
properly understood by perturbing the Einstein equations of general relativity a round the 



homogenous solution described in §4.1.31 - the details of this can be found in, e.g., iDodelson 



lOgi 

(|2003l 1. Here we go through a Newtonian approximation that elucidates much of the key 



physics governing cosmological perturbations, and then discuss the statistical measures used 
to analyze them. 



Linear growth of structure 

The equations governing the behavior of a non-relativistic, self-gravitating fluid are the 
continuity equation, Euler's equation, and the Poisson equation, respectively given by 

§^ + Vr-(pv) = (4.13) 

5V Vr-P 

+ (v.Vr)v = ^-Vr$ (4.14) 

at p 

Vl^ = AnGp (4.15) 

where p (r, t) is the density of the fluid, P (r, t) is the pressure, v (r, t) is the velocity of 
a fluid element, $(r, i) is the gravitational potential., and Vr denotes the gradient with 
respect to r. The spirit of first order cosmology is to start with a homogenous universe and 
then consider small perturbations to p, v, $,and P. We will use 

$o(r,t) = ^Po(r,t)r2 
Po(r,t) = 

for our homogenous solution, which satisfies equations (14.130 . ()4.14p . and ()4.15p if a(t) 
satisfies our zeroth order cosmology equations (j4.3p and (j4.8p . This corresponds to a universe 
of pressureless matter (such as dark matter) undergoing Hubble expansion vq = Hr. 
Now we solve for the behavior of small perturbations: 

■ po {r,t)+pi (r,t) 

= vo(r,t) + vi(r,t) (4.17) 
= $0 (r,t)+$i(r,t) 

where quantities with 1 subscripts are small enough that we can ignore all but the linear 
terms. Inserting these into the continuity, Euler, and Poisson equations, subtracting the 
homogenous solution, and dropping higher order terms of the perturbations gives 
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^ + (vo • Vr) p + poVr • VI + 3-pi = (4.18) 
ot a 

^ + (VO • Vr) VI + -VI = -Vr^i (4.19) 

ot a 

Vl^i = AttGpi. (4.20) 

Cosmology is often best understood in comoving coordinates defined by x = r/a{t). 
Comoving coordinates are defined such that they remain fixed for observers moving along 
with the Hubble expansion in a homogenous universe. The comoving peculiar velocity is 
given by u = vi/a(t), and we define the density contrast as 5 = pi/po- Re- writing the 
above equations in comoving coordinates gives 



5 = -V-u (4.21) 
U + -U = -^V$i (4.22) 
V^^i = 47rGaVo5, (4.23) 
where V = a (t) Vr is the gradient with respect to the comoving coordinates x, dots repre- 

1 /2 

sent d/dt = + vq • Vr, the time derivative for comoving observers, and Cg = {dP/dp) ' 
is the sound speed. 

By taking the time derivative of equation ()4.2ip and the divergence of equation (|4.22p 
we can derive a second order differential equation for 5: 

5 + 2-5- ^■kGpo5 = 0. (4.24) 
a 

This is the basic equation governing the growth of small density perturbations in a universe 
of pressureless matter, which is quite applicable to our universe since the bulk of large-scale 
structure formation in our universe occurs while dark matter is the dominant component. 
Equation (|4.24p has two linearly independent solutions: 

5g (t) = f ^da (growing mode) .^ ^5) 

(i) = f (decaying mode) 

where the growing mode increases with time and the decaying mode decreases. The general 
solution is a linear combination: 

5 (x, t) = A (x) 5g {t) + B (x) 5d it) , (4.26) 

where A and B are arbitrary functions of x. The decaying mode will eventually become 
negligible as the universe expands, and is thus not relevant to structure formation, so we 
focus on the growing mode. We can then write the density contrast as 

5(x,t) = 5(x,t,)^, (4.27) 

where the perturbations are specified at some initial time ti. Thus an initial perturbation 
field (5(x, tj) will retain its shape but increase in amplitude until 5^1 and our linear 
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approximation no longer holds. This is the basic picture of how small perturbations start 
growing into large structu res through gravitatio n al instabili t y. For a more detailed, properly 
relativistic treatment, see IPadmanabha^] (|l993l ): lOodelsonl (|2003l l. 



Statistics of random fields 

A general theory of how structure formed in our universe could never hope to predict the 
value of the density contrast 5 at a specific point x. Instead the goal is to predict its 
statistical properties, which are typically described in the language of random fields. A 
random field is an infinite dimensional random variable that has some probability to take 
on a given value at every point in space. The function 6{x,ti) in equation (14.27p . for 
example, is a realization of a random field. 

The statistical properties of a random field 6 are specified by a series of joint probability 
distributions Pn {8i, 62, ■ ■ ■ 6n) for its values 5i = (5 (xj) at n points xi,X2,...x„ for n = 
1, 2, 3, . . . . The most commonly used random field in cosmology is a zero-mean Gaussian 
random field, where the probability distribution is a multivariate Gaussian: 

Pn {5l,52, ...5n) dSid52 . . . d5n = ^, cxp if 8) d5idS2 . . . dSn, (4.28) 

Y^(27r)"det(M) 

where (5 is a column vector of the (5j's and M is a covariance matrix that defines the distri- 
bution. The expectation value of any function / of the Jj's is defined as 

if (61,62,... 6n))= J f {61,62,... 6n)Pn {61,62,... 6n)d6id62...d6n. (4.29) 

In typical cosmological applications, we try to infer the statistical properties of a ran- 
dom field, e.g. the density contrast, from observing the one realization of it provided by our 
universe. We do not have the luxury of observing multiple universes drawn from the ensem- 
ble defined by our model probability distribution, so we cannot really measure an ensemble 
average like equation (j4.28p specifies. To connect our statistical model with what we can 
actually observe, we have to make an assumption that the random fields are ergodic, which 
means that averaging over a sufficiently large volume of space is equivalent to averaging 
over a statistical ensemble, so we can use 

if {61,62,... 6n)) = lim / / ... / f{6{xi),6{x2),...6{xn))d^xid^X2...d^Xn 

(4.30) 

for a practical means of defining expectation values. 

The most commonly used statistical measures of random fields are moments of the 
probability distribution, called n-point functions. The 1-point function, i.e. the mean, is 
given by {6i), the 2-point function is given by {6i6j), and so forth. For example, a zero- 
mean Gaussian random field with P„ given by equation (j4.28p has {6i) = (hence the 
"zero-mean") , {6i6j) = Mij, and {6162 ■ . . 6n) = for n > 2. 

The 2-point function is the most basic statistical measure of inhomogeneities in a random 
field and thus plays an important role in first-order cosmology. It is known as the correlation 
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function ^ (xj,Xj) = {5i6j), and for a universe that is statistically homogenous and isotropic, 
the the correlation function will only depend on the distance between the two points Xj and 
Xj and not on their absolute position or orientation: ^ (xj,Xj) = ^ (|xj — Xj|). 

The 2-point function has some particularly convenient properties if we transform into 
Fourier space. Defining the Fourier transform as 



^ (k) = J e'^ ''6 (x) d^x, (4.31) 
the Fourier space version of the 2-point function is 

(^6* (ki) 6 (k,)) = I e*k,■x,-^k,.x.^ (^^^ ^ .) ^3^^^3^ (4 32) 

and if the correlation function depends only on the separation, i.e. (xj, Xj) = (xj — Xj), 
then equation (j4.32p reduces to 



5* (k,) 5 (k,)^ = (27r)^ 5d (ki - k,) P (k^) , (4.33) 

where 5_d is the Dirac delta function and P (k) , called the power spectrum, is the Fourier 
transform of the correlation function: 

P{k) = J e^'^-^'e (x) d^x. (4.34) 

Furthermore, if the correlation function depends only on the magnitude of the separation, 
i.e. ^ (xj,Xj) = ^ (|xj — Xjl), then the power spectrum will depend only on the magnitude 
A; = |k| of the wavevector k: P (k) = P (k). 

Equation (j4.33p illustrates why theorists prefer to work with the power spectrum rather 
than the correlation function: the Dirac delta function indicates that for a random field 
that is homogenous in real space, its Fourier modes will be uncorrelated for different values 
of k. In other words, its covariance matrix is diagonal: if 6 (x) is a homogenous, zero- 
mean Gaussian random field described by the probability distribution in equation (j4.28p . 
then 6 (k) is a Gaussian random field as well, with the additional feature that the matrix 
M in equation (j4.28p is a diagonal matrix. This makes the pow er spectruin muc h more 
convenient for many theoretical calculations. For more details see Hamilton ( 19981 ) for an 
excellent discussion of correlation functions and power spectra. 

The primary goal of first order cosmology is to predict the power spectrum P {k, t) of 
the density contrast 6 (x, t) of matter in our universe as a function of cosmic time. This can 
be modeled as 

p{k,t)=(j^yT'{k)P{k,u), 

where 6g (t) is (an appropriately relativistic version of) the growing mode in equation (j4.25p 
and T (k), known as the transfer function, accounts for a bunch of physics we have ignored 
in the simple treatment here - interactions between matter and radiation, the behavior of 
modes with wavelengths larger than the horizon (the distance a pho t on can travel through 
the universe since time t = 0) and so forth. See Padmanabhan ( 19931 ) for a detailed 
discussion of the transfer function. 

The function P (k, ti) is the spectrum of primordial fluctuations that seeds the growth of 
structure. Theories of inflation predict that this initial spectrum was produced by quantum 
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fluctuations that were stretched to macroscopic size by a period of exponential expansion 
in the very early universe. This process typically predicts that the initial fluctuations can 
be described by a Gaussian random field with a power-law power spectrum: 



P (A;, tj) — Agk " , 

where Ag and are the basic first order cosmological parameters that are included along 
with the zeroth order parameters Hq, w, and various fi's in parameter estimation analyses 
such as the one discussed in §4.31 Most models of infiation predict that to be slightly 
less than 1. 

Theoretical predictions for P (k, t) for matter in our universe given a set of cosmological 
parameters can be computed exactly i n the linear regime and is ty pically done using the 



handy computer program CMBFAST (ISeljak and Zaldarriagalll996l ). With the first order 



approximation that 5^1, each Fourier mode evolves independently, so this makes the 
theoretical calculations much easier in Fourier space than in real space. However, there are 
two issues with comparing this predicted power spectrum with what we observe: first, for 
perturbations with (5 ~ 1 the linear theory no longer applies. Secondly, our observations rely 
on using galaxies to trace the overall matter distribution. The issue of nonlinear evolution 
is discussed in the next section, and the issue of galaxies as tracers is explored in depth 
in chapter [6l Here we merely note the usual assumption made when measuring the power 
spectrum from galaxy surveys: that the power spectrum of the galaxies has the same shape 
as that of the overall matter, but scaled by a normalization factor: 

^galaxy (k) = 6'Pmattcr (k) , (4.35) 

where the parameter b is called the bias. On scales larger than the size of the largest grav- 
itationally bound structures in our universe {k < O.Og/i/Mpc, or A = 27r/A; > TO/i^^Mpc), 
the simplifying assu mptions of linear evo lution and simple galaxy bias have been found to 
be quite reasonable ( Tegmark et al. 20061 ). 



4.1.5 Higher order cosmology 

The zeroth and first order cosmology discussed in the previous sections can take us a long 
way, but a full theory of structure formation in our universe requires a nonlinear treatment 
of differential equations such as equations (|4.13p . (|4.14p . and (|4.15p that can be applied to 
systems like galaxies and galaxy clusters that have 6^1. Unfortunately, this turns out to 
be quite challenging - the nonlinear versions of the equations couple together the evolution 
of different Fourier modes, so most of the techniques used in first order cosmology cannot 
be applied and exact solutions do not exist. 

This difficult and unsolved problem has been tackled from a number of different angles: 
analytical models that apply under various simplified conditions, numerical simulations 
that attempt to include all of the relevant physics, empirical fitting formulas tuned to 
match simulations or observations, and various hybrid techniques that endeavor to combine 
analytical and numerical results to build up a complete picture. Numerical simulations 
are arguably the only way to obtain a complete solution, and much w ork has been done on 
both A^-body simulations that model purely gravitational systems fe.g. Springel et al. 20051 ) 



and hydrodynamical simulations that include gas physics as well (see iDolag et al.ll2008l for a 
recent review). However, such simulations are al ways subject to a variety of approxiination s 
and the results can be challenging to interpret ( Tasitsiomi et al.l 2004; Romeo et al. 20081 ). 
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so analytical and semi-analytical methods continue to play an important role. Here we 
outline perhaps the most widely-used analytical model for nonlinear structure evolution 
and discuss some of its extensions, and then we give an overview of a hybrid approach 
known as the halo model that has had much success in recent years. 



Spherical collapse and Press-Schechter theory 

The analytical model we will discuss here starts with a simple approximation for the collapse 
of a spherical overdensity. Consider a universe that contains only matter and has a density 
equal to the critical density pcr such that 0,^ = 1, except in a spherical region of comoving 
radius Rq enclosing mass M = 47rr2^/3crO-Ro/3 in which the density is slightly higher than 
the critical density, with > 1. Equation (j4.12p can be solved for the = 1 background 
universe to give 

/ 2 N 2/3 -j^ 
at (t) = i^-Hotj , pb {t) = (4.36) 

where the b subscript denotes the background universe. The perturbed region also evolves 
according to Equation ()4.12p . but with Q'^ > 1. A parametric solution is given by 

«p(^) = ^n!?^/l-cosr?) ^^^^^ 
^W = 25o (nJ->^ (r?-sin7?) 

where the p subscript denotes the perturbation. The perturbation will expand until r/ = vr, 
and then start to contract until r] = 2tt, at which point it will be fully collapsed. 

In reality, the density will not become infinite - if the matter is baryonic, the collapse will 
eventually be stopped by gas pressure, and if it is dark matter, the particles will whiz past 
each other and get slowed down by gravitational forces. Either way, the system is roughly 
expected to reach an equilibrium state where its kinetic and potential energy should satisfy 
the virial theorem: U + 2K = 0. The virialized radius can be found by considering the 
energy at the turnaround point r/ = vr is purely potential. Since the gravitational potential 
energy of a sphere is proportional to R^^, energy conservation and the virial theorem imply 
that the virialized radius is half the radius at turnaround, so the collapsed density will be 

/OcoUapsed = 8pp {l] = Vt) 

= ISvrVb iv = 2vr) . (4.38) 

After the collapse, the virialized object will retain this physical density as the background 
universe continues to expand. 

The density contrast between the perturbation and the background universe can be 
found from equations (|4.36p and (j4.37p to be 

^^Pp-Pb _ 9{r]-smr]f 



Pb 2(l-cosr/)'^ 

~ ^T]"^ fovr] -^l. (4.39) 

How does this compare to the first order evolution of 6 derived in §4.1.41 .'' Equation 14.251 
for a flat, matter-dominated universe gives 6g cx a, so the linear theory prediction gives 
(5iin oc i^/^. Choosing the proportionality constant to match equation 14.391 for early times. 
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lin = 7;7T[6(r/-sm7?)]'/^ (4.40) 



we have 

3_ 
20 

Equations ()4.39p and (|4.40p give us an important tool: a means of connecting an exact solu- 
tion with linear theory predictions. We have done these calculations fo r a the simple rnodel o f 
a = 1 universe, but the same basic ideas apply for other cases - see PadmanabhanI ( 19931 ) 



for more details. At the time of the collapse, we have (5iin (r? = 27r) = 3 (127r)^/^ /20 « 1.686. 
Based on this, we can construct a recipe for modeling nonlinear evolution: evolve 5 using 
linear perturbation theory until it reaches a critical value of 6c = 1.686 in a given region, 
and then substitute in a spherically-collapsed, virialized object at that location. 

Press and Schechter ( 1974j ) applied this idea to describe structure formation in our uni- 
verse and produced a canonical formula for the mass function n (M) giving the comoving 
number density of objects of mass M based on the following argument: consider a uni- 
verse with mean comoving mass density p that has a matter distribution described by the 
linearly-evolved power spectrum P^^ (k). To determine the properties of the population of 
mass M objects, smooth the linearly-evolved density field with a spherical top-hat filter 
w (r, R) = 3/{4ttR^)Q{R - r) with radius R = (3M/ (47rp))^/^ The smoothed density field 
is 

§;^j (x) = j w{\x- x'l , R) 6 (x ) d^x, (4.41) 
and the variance of the smoothed field is 

{5m i^r) ^ a^M) = I f ^%^|^(^,ii)P, (4.42) 

where w {k, R) is the Fourier transform of the top hat filter: 

3 

w{k,R) = [sin (kR) - kR cos (kR)] . (4.43) 

Now consider any region that has Sm (x) > 6c to be a virialized object that has mass > M. 
By repeating this for all smoothing scales, we can build up a picture of the mass distribution 
of objects in our universe. 

Assuming that the density contrast 6 can be modeled as a Gaussian random field, the 
fraction of the mass in our universe contained in virialized objects with masses > M is given 
in the Press-Schechter model by 

1 f°° / S"^ \ 

F{>M) = - / exp --^^ d,5M 



^erfc (.[^] , (4.44) 



2 VV 2^ 

where v = 6c /c'^ i^)- This formula has a curious feature: it says that the maximum 
fraction of mass that can be contained in virialized objects is 1/2, not 1 - this is because we 
have ignored underdense regions that have 6 < 0. However, mass in under dense regions is 
likely to accrete onto the virialized objects eventually. To account for this omission, Press 
and Schechter simply multiply equation (j4.44p by a "fudge factor" of 2 and find they get a 
sensible final answer. 

We can find the mass function by differentiating equation (j4.44p and multiplying by the 
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extra factor of 2 and by the number density p/M that would be expected if all of the mass 
in a universe with average mass density p were in mass M objects: 



n (M) dM 



MdM 




dM 



(4.45) 



This is the famous Press-Schechter mass function. Although it is based on a rather heuristic 
argument, it has proven to be a rem arkably valuable tool - it agrees well with numerical 
simulations (e.g. White et al. 19931 ) and has provided a launching point for further theo- 
retical work, p articular l y the extended Press-Schechter theory based on the excursion set 
formalism (s ee IZentnerl (120071) for a recent review) . Ear ly applications of the excursion set 



formalism (jPeacock and Heavenslll990l : lBond et al.lll99ll ) provided a more formal derivation 
of equation ()4.45p that explains the factor of 2 and resolves other issues as well. 

An upda ted version of the ma ss function t hat provides a bette r fit to the latest numerical 
simulations ( Springel et al. 20051 ) is given by Sheth and Tor men ( 19991 ): 



n (M) dM 




exp 



qv 



d{\nu) 
d (In M) 



dM, 



(4.46) 



with q « 0.707, p « 0.3, and A{p) = {I + 2~pT {1/2 - p) / ^) ^ k, 0.3222. Note that 
i f p = and q = 1 then This mass function can be derived using excursion set theory 
(jSheth et al.ll200ll ) yy considering ellipsoidal rather than spherical collapse. 

The extended Press-Schechter theory also provides a wa y to make theoretica l predi ctions 



see 



Martinez and Saar (j2002l ) for a 



for the clustering of halos an d their merger his tories , , 

good summary. For example, Bond et al. ( 199ll ) showed that the fraction of mass of a large 
object of mass Mq with present linear density contrast 5q and variance ctq that is made up 
of merged- in smaller objects of mass Mi, present linear density contrast 6i and variance ai 
is given by 



/(cJi,5i|cJo,5o)— f dMi 
dM-[ 



So 



2tt 



exp 



{Si - Sof 
2 {al - ai) 



^dMi 



dM 



(4.47) 



Mo and White ( 19961 ) use equation (j4.47p to calculate the bias factor bM for halos of 



mass M using the extended Press-Schechter theory, where is defined such that the power 
spectrum Pm {k) of virialized objects of mass M is f'|,/-Piin {k), with Piin {k) being the linear 
theory matter power spectrum. They find that 



1 + 



1 



(4.48) 



The generalized version of this equation using the Sheth- Tormen mass function (equa- 
tion ()4.46p ) was derived by Sheth and Tormen ( 19991 ): 



bM = l + 



qu — 1 



+ 



2p/Sc 



(4.49) 



Sc ' l + {qi^r' 

These and other results from extended Press-Schechter theory are key ingredients in assem- 
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bling an analytical description of nonlinear structure formation. 



Halo model 



One particular model that has proven useful for describing both nonlinea.r clustering and 
galaxy bias is the known as the halo m odel, first formalized by lSelja IScoccimarro et al 

(|200lh . rSee ICoorav and Shethlliooi ibr a comprehensive review.) The basic idea is this: 
at late times, the matter distribution in our universe can be modeled using only virialized 
objects made of dark matter that are called halos, and the galaxy distribution can be mod- 
eled by an appropriate scheme for placing galaxies in these dark matter halos. The key 
assumption of the model is that the average properties of a halo - its density profile and its 
galaxy content - only depend on the halo mass. 

Using this assumption, we can construct an analytical form for the nonlinear power spec- 
trum by assembling components from first order cosmology, the extended Press-Schechter 
theory discussed above, and A^-body simulation results. The halo model power spectrum is 



P (k) = P^^ (k) + P^^ (k) , (4.50) 

where P^^ (k) is the 1-halo term arising from correlations between mass within the same 
halo and P^^ (k) is the 2-halo term arising from correlations between mass in different halos. 

To calculate these terms we need a descri ption of how rnass is distributed within a halo. 
This is usually modeled by an NFW profile (INavarro et al.lll996l '). which has been found to 
be an excellent fit to dark matter halos that form in A^-body simulations. The normalized 
profile u (r|M) = p{r\M) /M is given by 



u{r\M) 



47vr 



vir / cr 

vir 



1 + 



(4.51) 



where f = 1 / [In (1 + c) — c/ (1 + c)] and the concentr ation parameter c is a function of the 
mass, typically approximated as (iBullock et al.ll200lh 



c(M, z) 



9 



M 



l + z \M^(z 



-0.13 



(4.52) 



with the characteristic mass scale is defined as the point where v = 1: 5l = a'^ (M*). 
The virial radius r^\^ is defined by the mass: M = A/347rr^jj./3, where A is the typical 
overdensity of a virialized structure compared to the background density - we found in 
equation (j4.38p that A = ISvr^ ~ 178 for a Q.m = 1- Equation ()4.5ip is normalized such 
that the total mass within the virial radius is equal to M. 
The 1- and 2-halo terms are then given by 



pi'^(A;) = ^ j dMn{M)\u{k\M)[ 



(4.53) 



P^f" (k) = ^ j dMin{Mi)u{k\Mi) j dM2n{M2)u{k\M2)P^^{k\Mi,M2), (4.54) 

where u {k\M) is the Fourier transform of the NFW profile in equation (j4.5ip and n (M) is 
the mass function given by equation (14461) . P^^ (A;|Mi, M2) is the halo-halo power spectrum. 
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most simply modeled by 



P''^ {k\M^,M2)=bMMhPlin{k), 



(4.55) 



where 6m is bias for halos of mass M gi ven in equation (I4.49D. although more sophisticated 



models have been used in recent work (jZehavi et al.ll2005l : [Tinker et al.ll2005l ) 



In addition to modeling the dark matter distribution, the halo model is also a useful 
tool for modeling galaxy bias. The galaxy distribution then be modeled by first deter- 
mining the halo distribution and then populating the halos with galaxies. The second 
step requires a model for P {N\M, L,C, . . .), the probability that a halo of mass M will 
host galaxies with luminosity L, color C, and any other galaxy properties the modelers 
wish to i nclude. Multiple parameterizations for this probability distribution have been de- 



20051 : lYang et al.l 



200: 



veloped (iPeacock and SmithI l2000l: iBerhnd aiid Weinberg! 1200 



van den Bosch et al. 



ty Qistribution nave been de- 
3; Sefusatti and Scoccimarrol 



20031 ) . The algorithm used for placing the 



galaxies in halos encodes the physics of galaxy for mation and can be used to compare semi- 
analytical galaxy formation models (|Baughll20od 'l with observational results such as those 
discussed in chapter ([6]). This procedure yields a model of galaxy bias that is more closely 
linked with the actual physical processes than the simple assumption of equation (j4.35p . 



The halo model has been successfully used to mode ^ 



tions , such as galaxy correlation functio ns (IZehavi et al 



seve ral different typ es of observa- 
2OO5I I. void statistics (iTinker et al 



20071 ) ■ and counts of galaxies in groups (lYang et al.ll2008l ). However, as we will see in §4. Sl it 
has yet to be incorporated into full cosmological parameter analyses - this could potentially 
be a promising direction for future work. 



4.1.6 What does this teach us about particle physics? 

What can we hope to learn about physics on at the sub-atomic level by studying our universe 
on such a grand scale? Quite a lot, as it turns out! One basic tenet of particle physics is 
that higher energies can probe shorter length scales. In order to examine extremely small 
length scales - the quarks within a proton, for example ~ particles are collided with each 
other at extremely high energies. The current generation of particle accelerators can achieve 
energies of 2 TeV, and the Large Hadron Collider, set to turn on in 2008, will push this up 
to 14 TeV ( Feng 20081 ). However, the extremely hot, dense conditions in the first fraction 



of a second after the Big Bang provide a natural particle accelerator - a split-second after 
the Big Bang, the average thermal energy of particles could have been as high as 10^^ TeV. 
Thus by studying cosmology, particle physicists can learn about phenomena that could not 
be created in a terrestrial experiment. 

Of course, the problem with this plan is that after 14 billion years, much of the action has 
faded away - unstable particles created early on have long since decayed and the present-day 
density of any remaining stable particles has been so diluted by the expansion that they 
are no longer interacting with each other. Particle physicists who endeavor to use our early 
universe as a laboratory must find ways to predict detectable signatures left behind by the 
processes they hope to investigate. Particle physicists also work to develop models that can 
explain the bizarre features of the standard cosmological model described in ^4.1.21 This is 
especially interesting because the ideas of dark matter and dark energy cannot be explained 
within the framework of the standard model of particle physics and thus hint at physics 
beyond our current understanding. 

From cosmological and other astronomical observations, we know that some form of 
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non-baryonic dark matter that does not couple strongly to photons makes up most of the 
mass in our universe. Theorized particles that have properties similar to dark matter are 
actually well- motivated in a particle physics context. One class of proposed extensions 
to the standard model of particle physics, known as supersymmetry, postul ate that ever y 
particle has a supersymmetric partner particle whose spin differs by 1/2 ( Peskin 20081 ) . 
The supersymmetric partners of known particles are expected to have such high masses 
that we would not have been produced by the current generation of particle accelerators. 
Supersymmetry was proposed to explain a number of seemingly "unnatural" features of 
the standard model, such as why the masses of standard model particles are so much 
smaller than the natural energy scale of the theory. It also provides a candidate particle for 
dark matter: the lightest supersymmetric particle. At the extremely high temperatures and 
densities present in the early universe, an abundance of supersymmetric particles might have 
been created. Most of these will be unstable and decay quickly into lighter supersymmetric 
particles, but the lightest supersymmetric particle is predicted to be stable. If this is the 
case, it would be present in our universe today and act like dark matter - that is, like 
a weakly interacting massive particle, a.k.a. a WIMP. This is currently regarded as the 
most natural explanation for dark matter, but there are a number of other particle physics 
theories that can generate p a,rticles with da rk-matter-like properties as well, such as axions, 
solitions, or WIMPZILLAS (|Gondololl20o3 ). 

One might ask why we need to resort to such exotic-sounding models to explain dark 
matter - after all, we know that neutrinos must have some mass, as discussed in ^ so 
could the dark matter simply be neutrinos? The answer comes from first order cosmology, 
considering the nature of the perturbations in the density field. The three types of neutrinos 
in the standard model, while we know they must have some mass, are still light enough 
that their average velocity is quite high at the time when structure is forming, so they are 
classified as hot dark matter. The main effect of hot dark matter is to wash out the formation 
of structure on small scales - heuristically, neutrinos will not collapse into a structure if 
their average velocity exceeds the escape velocity of the a given overdense clump. For small 
(i.e., galaxy-sized) clumps, the escape velocity is sm all enough that we would not e xpect 
neutrinos to collapse into structures of this size. See iLesgourgues and Pastoil (|2006l ) for a 
detailed explanation of this effect. 

Such a picture is not consistent with the structure we see today, which implies that the 
dark matter must be cold - that is, massive enough to have very low average velocity at 
the time structure starts to form. In fact, measurements of the power spectrum can be 
used to set a strong upper limit on the sum of the masses of the three types of neutrinos 
- we discuss an example of this in §4.31 In fact, cosmology currently p lace stronger upper 
limit s on the neutrino mass than current particle physics experiments (|Elgar0v and Lahav 
2006 : Hannestad et al. 20081 ). and the next generation of cosmological measurements may 



be sensitive enough to detect masses as small as the mass differences seen by Super-K. 

Another cosmological puzzle particle physicists have sought to tackle is that of dark 
energy. What kind of substance could lead to the bizarre repulsive gravitational effects that 
are driving the current acceleration of our universe's expansion? The most basic possibility 
from the point of view of quantum field theory is that there is an inherent energy density 
associated with empty space called vacuum energy. The vacuum energy density would 
remain constant as the universe expands, which is consistent with the observed properties of 
dark energy. Theoretical predictions of what the value of the vacuum energy density should 
be are rather problematic - a basic quantum field theory calculation gives a divergence, 
predicting infinite vacuum energy density. Imposing a natural cutoff to this divergence 
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gives an estimate of ~10^^'^ erg/cm^ for the vacuum energy density. Prior to the discovery 
of dark energy, particle theorists devoted much effort to searching for natural ways this 
vacuum energy could cancel out to be exactly zero. Finding that there is indeed something 
that looks like vacuum energy in our universe made the problem even more difficult for 
the theorists, since the cosmological observations indicate that the density of dark energy 
is ~10~^'' erg/cm^, 120 orders of magnitude less than the estimated value! A theory where 
the vacuum energy is almost entirely canceled out except for a tiny fraction would be much 
more bizarre than perfect cancellation. Supersymmetry provides some help here, but not 
enough - it can reduce the discrepancy to "just" 60 orders of magnitude. D ark energy has 



prove n to be a fascinating challenge from a particle theory perspective (see lFrieman et al 



20081 for a review of possible models) and is particularly interesting from the point of view 
of connecting particle physics and astronomy because as far as w e know, the only way to 
measure its properties is through cosmological observations (see lAlbrecht et al.l l2006l for 
more details). 



4.2 Tools of cosmology 

4.2.1 Combining cosmological probes 

The great strength of modern observational cosmology is that such a wide variety of differ- 
ent experimental measurements point to the same theoretical model. Some measurements, 
such as using supernovae as standard candles to trace the expansion history of our uni- 
verse, target the zeroth order cosmology discussed in ^4.1.31 - that is, they are effectively 
probing the function a{t). Other measurements, such as studying the tiny fluctuations in 
the cosmic microwave backround or the clustering patterns in the distribution of galaxies, 
are fundamentally probes of the first order cosmology discussed in §4.1.4t they measure 
variations on the power spectrum P{k). Still others, such as measuring the abundances 
of elements formed by fusion in the hot, dense plasma of the early universe, target one 
particular cosmological parameter - in this case, the ratio of baryon density to photon 
density. 

Combining all of these complementary probes is an extremely powerful way to study 
our universe. First of all, since different probes are sensitive to different combinations of 
cosmological parameters, combining two or more different measurements allows us to break 
the degeneracies between the parameters we want to measure. Furthermore, because each 
technique is subject to different systematic effects, we can use them as consistency checks. 
If two experiments disagree, it would be a sign that either there is some systematic effect 
that has not been accounted for or that our theoretical model is flawed. The fact that there 
is broad agreement over so many different experiments lends much stronger support to the 
standard cosmological model than any single experiment could ever do. In the following 
sections we describe two key probes - the cosmic microwave background and galaxy surveys 
- in some detail, and give a brief overview of other complementary techniques. 

4.2.2 Galaxy surveys 

The work in chapters [5] and [6] focuses on issues with using such surveys as cosmological 
tools, so we describe this particular cosmological probe in some detail here. The basic idea 
behind a galaxy survey is to map out the locations of galaxies over a large volume and use 
them as point tracers of the overall matter distribution in our universe. This provides a 
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way to observe the statistical properties of the large-scale structure predicted from the first 
order and higher order cosmological models discussed in §4. 1.41 and §4.1.51 

Galaxy surveys come in two basic flavors: angular surveys and redshift surveys. Angular 
surveys indentify galaxies over a large area of the sky, and thus only see the two-dimensional 
projection of the galaxy distribution. Redshift surveys, on the other hand, also measure 
spectra to obtain the redshift of each galaxy and thus its distance from us according to 
Bubble's law (equation (14. 4p ). which produces a three-dimensional map of galaxy loca- 
tions. Recently a hybrid approach between these two types of surveys has become popular: 
photometric redshift surveys, whereby the redshift is estimated based on the photometric 
information from multiple wavelength bands of an angular survey. This provides rough 
measure of the distance that eliminates the need for time-consuming spectroscopy and the 
accuracy of the redshift estimation tech niques is continuing to improve (ICollister et al.ll2007l : 



Ovaizu et al.ll2008l : iBanerii et al.ll2007l l. 



Furthermore, surveys can either be wide and shallow - covering a large area of the sky 
but only extending a modest distance from us, or narrow and deep - focusing on a small 
patch of sky for an extended time to find very faint and therefore extremely distant objects. 
Typical surveys strike a balance between these two choices that is optimal for a particular 
set of science goals. 



Early observations 



Astronomers have been interested in studying th e distribution of galaxies since the time of 
Hubble (see, e.g., Hubblel 1934 : Carpenter 1938). The early angular catalogs studied were 
the Palomar ( Abell 1959 ) and the Lick ( Doughty et al. 1974 ) sky surveys, and redshifts of 



over 8 00 galaxies from various spectroscopic observations were collected in iHumason et al 
(|l95(il ). These early studies showed evidence for dumpine ss in the overall distribution 



Zwickvl I1952I: Ide Vaucouleurs 



although its exac t nature was difficult to discern (see, e.g 

19581: lAbelll Il96ll ) These observations inspired the theoretical work of iNeyman and Scotd 
(1952) which formed the seeds of the present-day halo model discussed in ^4.1.51 as well as 



the f irst effort to compare observations with a simulated galaxy distribution (IScott et al 
I954I ). 

The first true redshift survey was the Center for Astrophysics Redshift Survey (CfA; 



Huchra et al.lll983l). whic h started in 1978 - prior to this only ~ 2000 galaxies had measured 



redshifts (|da Costal [l999l ) . This survey measured redshifts of the 2401 galaxies with an 
apparent Z-band magnitude brighter than 14.5 m^. A slice through the CfA survey is 
shown in Fig. 14-21 This survey revealed that galaxies are distributed in a complex web of 
filaments, sheets, and voids. 

Over the following decade a numb er of other surveys w ere undertaken, including the 
Southern Sky Redshift Sur vey (SSRSl; da Costa et al. 199ll ) , surveys done by the Infrared 



Astronomy Satellite (IRAS; Strauss et al. 19921 ). the Automatic Plate Measurement angular 



survey (APM; Maddo x et al. 1990a,b, 1996), and th e Las Campanas Redshift Survey (LCRS; 



Shectman et al. 19961 ) 



Peacock and Dodda (Il994l ) compiled power spectrum measurements 
from this generation of surveys and fit them with several cosmological models - they found 
that the data favored models with low matter density: ^Im < 1, in agreement wit h present 



day st a ndard rnqdel. Further details on the history of this field can be found in Ida Costa 
(jl999l ): lBiviM^ (|2000l ). 
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Figur e 4-2 A slice through the CfA redshift survey (jHuchra et al 



1^) done at the Smithsonian Astrophysical Observatory. Figure from 



http : //cf a-www . harvard . edu/ "^uchra/zcat/ 



Modern surveys 



The current state of the art in galaxy survey cosmology is exemplified by two sur veys: the 
Two Degree Field Galaxy Reds hift Survey (2dF GRS; [Colless_etal| [200l|, l2003l ^ and the 
Sloan Digital Sky Survey fSDSS: lYork et al.ll2nnd l. 

2dFGRS is a redshift survey that measured the redshifts of galaxies observed by the 
APM angular survey. Target galaxies with an apparent blue-band magnitude brighter than 
bj < 19.5 were selected from APM's 5° x 5° photographic plates and spectra were taken 
using the Two Degree Field system at the Anglo- Australian Telescope in New South Wales, 
Australia. This is a 4 meter telescope with a 2° field of view and a set of spectrographs 
which accept 400 optical fibers such that spectra of 400 objects in the field of view can 
be taken simultaneously. The redshift survey was completed in 2002 and took 272 nights 
over 5 years. It measured redshifts for 221,414 galaxies with a mean redshift of 2; = 0.11 
and covered an area of about 1500 square degrees. A slice through the galaxy distribution 
measured by 2dFGRS is shown in Fig. 14-31 

SDSS, our source of data for the work in chapter [6l is the largest galaxy redshift survey 
to date, and is continuing to gather more data. It is being conducted using a dedicated 2.5 
meter telesc ope at the Apache Point Observatory in New Mexico that has a 3° field of view 



and ~ 1.4" (|Gunn et al.l l200(il ^. The final SDSS data set will cover nearly 1/4 of the sky 



and have spectra of nearly 1 million galaxies and about 105,000 quasars (a.k.a. AGNs; see 
§2.3.2p . The most recent data release, DR6, covers 17% of the sky and has 790,860 galaxies 
and 103,647 quasars. 

The camera used for the photometric survey is a 5 x 6 array of CCDs arranged such that 
each of the 5 rows ob serves through a different filter, ranging from near-ulraviolet to near- 
infrared wavelengths ( Fukugita et al. 19961 : Gunn et al. 19981 ). The survey method used is 
known as drift-scanning: as the sky moves overhead due to the rotation of the Earth, points 
on the sky pass through each filter in succession. This produces a long strip of 6 scanlines 



75 




across the sky with photometric data from ah 5 filters simultaneously. It takes over 8 hours 
to scan a strip 130° long. 

Based on this photometric information, targets are selected for the spectroscopic survey 
to obtain redshifts. There are several types of targets selected, including "main" galaxies, 
the luminous red galaxies (LRGs), quasars, and stars within our galaxy. The main galaxy 
sample includes all galaxies with an r-band apparent magnitude brighter than r < 17.77 - 
a typical galaxy at a redshift of z = 0.1 has r = 17.5. The brightest galaxies in the main 
sample can b e seen out to a redshift of z ^ 0.3. This is the sample used in chapter [6l The 
LRG sample ( Eisenstein et al. 200ll ) is a sparser sample extending to higher redshifts that 
contains only galaxies of a particular type: very luminous early-type elliptical galaxies. The 
LRGs is roughly volume-limited and extends to a redshift of z ~ 0.5. The properties of the 
LRG sample make it especially useful for certain types of cosmological measurements, such 
as the study discussed in ^4.31 

SDSS spectroscopy is do ne using a pair of int egral field spectrographs that can take up to 
640 spectra simultaneously (|Uomoto et al.ll2004l ). This is done using 3° diameter aluminum 
plates drilled with holes at the locations of the targets and plugged with optical fiber cables 
which are fed into the spectrographs. Typically six to nine spectroscopic fields can be 
observed on a given night. Once a spectrum has been measured for a galaxy, its redshift 
can be determined by the positions of known spectral lines, and thus we know its three- 
dimensional location within the survey volume. A slice of the SDSS galaxy distribution is 
shown in Fig. 14-41 

Both 2dFGRS and SDSS have been extraordinarily valuable resources for the general 
scientific community, in large part because the data is made publicly available. The un- 
precedented volume mapped by these two surveys have revolutionized studies of large-scale 



76 



-24<M< 
-23<M< 
-22<M< 
-21<M< 
-20<M< 
-19<M< 
-18<M< 
-17<M< 



100 h"^Mpc 



Figure 4-4 A slice through SDSS (|York et alJbood ). The galaxies are colored according to 
their absolute r-band magnitude , using the same luminosity bins as in chapter [6l 



77 



structure in our universe in many ways. One worth noting here is that the galaxy samples 
are now large enough that they can be reasonably divided up into subsets for more detailed 
analysis, as is done in chapter [6l 

Other current surveys such as the 2-Micron All Sky Survey (2MASS: lHuchra et al.ll2005,| 



and the second phase of the Deep Extragalactic Evolutionary Probe (DEEP2; iDavis et al 
20031 ) the have led to quite interesting cosmological results as well - we have focused on 



2dFGRS and SDSS here because their balance on the wide and shallow versus narrow and 
deep spectrum is ideal for large-scale structure studies. However, the full-sky coverage of 
2MASS provides of our n earby universe has been quite fruitful for matching local velocity 
flows to the mass density (ICrook et al.ll2007l ') and the high redshift data f rom DEEP2 allo ws 
us to measure the evolution of galaxy properties and clustering as well ( Coil et al. 20081 ). 

There are a number of upcoming surveys planned that will take this field another giant 
step fo rward. Starting in the nex t year or so, the Baryon Oscillation Spectroscopic Survey 
fBOSS lSDSS Collaboration l2007l ) will use the SDSS telescope with updated technology 
to tar get the baryon oscillation signal describ ed in §4.2.41 and the Dark Energy Survey 
(DES; iDark Energv Survev CollaboratiOT]l2005l ) will perform a photometric redshift survey 
that will contain about 300 million galaxies. Further in th e future are the Panoramic 
Survey Telescope and Rapid Resp onse Systein (Pan-STARRS; Kaiser 20041 ) and the Large 
Synoptic Survey Telescope (LSST; Tyson 20021 : Stubbs et al. 2004), which will operate with 
unprecedented speed, surveying the entire available sky several times each month, and the 
Wide-Field Multi-Object Spectrograph (WFMOS; lYamamoto et al.l boOfil ). which will be 
able to take spectra of 20,000 galaxies every night. 

On the theoretical side, it is important to be prepared for these upcoming surveys by 
understanding systematic effects such as nonlinear structure growth and galaxy bias so that 
we can extract as much cosmological information from them as possible. Galaxy surveys 
will continue to be an important probe as cosmology moves into a new era. 



4.2.3 Cosmic microwave background 

Galaxy surveys map out the present-day or relatively low redshift universe, so they are par- 
ticularly well-complemented by probes of our universe's early stages. Perhaps the richest 
source of information we have about the early universe comes from the leftover radiation 
from our universe's hot and dense early phase, known as the cosmic microwave background 
(CMB). Measurements of the CMB are complementary to galaxy surveys that aim to map 
out our universe today because the CMB is essentially a snapshot of the early universe. 
Combining information from galaxy surveys and the CMB provides a way to break degen- 
eracies between parameters and tighten the constraints on our theoretical model. 

The basic physics behind this early universe snapshot is as follows: the early universe was 
filled with a plasma of protons and electrons coupled to photons in thermal equilibrium. 
The electrons were tightly coupled to the photons through Thomson scattering, and the 
protons are tightly coupled to the electrons through electromagnetic attraction, so these 
three components behave as a single fluid. 

The dark matter, however, does not interact with baryons or light, so at this point 
it has already started to clump together under the influence of gravity around the seed 
fluctuations generated by inflation. The mass of the protons drags the baryon-photon 
fluid into the gravitational potential wells formed by the clumping dark matter, but the 
radiation pressure from the photons resists this infall. This process generates oscillations 
which propagate as sound waves through the baryon-photon plasma. 
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Figure 4-5 Temperature anisotropies in the CMB measured by the NASA/WMAP Sci- 
ence Team, shown iii Hammer- Aitoff projection in Galactic coordinates. Reproduced from 
Hinshaw et al. ( 20081 ): used with permission. 



As our universe continues to expand, it also cools, and when our universe is about 
400,000 years old (0.003% of its current age), it became cool enough for the protons and 
electrons to combine into neutral hydrogen. At this point, the baryons and photons decouple 
from each other. The baryons are now free to collapse into the potential wells formed by 
the dark matter, and begin forming into stars and galaxies. The photons, on the other 
hand, free-stream through our universe and undergo (nearly) no further interactions from 
the time of decoupling until the present. 

Today these CMB photons permeate our universe, and they have been redshifted to 
microwave wavelengths as our universe has expanded. Because they have traveled largely 
unimpeded since they decoupled, they provide a snapshot of the early universe at the 
moment of decoupling. This leftover glow from the hot and dense early universe has been 
observed to be essenti ally a perfect blackb ody with a uniform temperature of 2.725K in all 
directions on the sky ([Mather et al.lll994l ). This fact alone is extremely solid evidence for 
the basic picture of the standard cosmological model described in §4.1.2[ 

However, the true power of CMB observations lies in tiny deviations from this unifor- 
mity: the CMB anisotropies. The sound waves propagating through the plasma of the early 
universe create slight inhomogeneities in the density, which we can see as tiny variations 
in the CMB temperature at different points on the sky. Figure 14-51 shows the latest mea- 
surement of these anisotropi es from the 5-year d ata release of the Wilkinson Microwave 
Anisotropy Probe(WMAP5; Hinshaw et al. 20081 ) . Since the fluctuations at the time of 
decoupling are very small (1 part in 10^), the first order cosmology described in §4.1.41 can 
provide extremely accurate predictions for the statistical properties of these anisotropies, 
namely their power spectrum. 
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Figure 4-6 CMB angular power spectrum measured by the NASA/WMAP Science Team 
from the 5 year data release. This plot shows the average m agnitude of the tempe rature 
deviations vs. their angular size on the sky. Reproduced from lHinshaw et al.l ( 20081 ): used 
with permission. 



To predict the power spectrum, we need only two basic physical assumptions: that we 
can apply general relativity to cosmological length scales and that our understanding of 
plasma physics is accurate under the conditions at the time of decoupling (a density of 
about 300 particles per cubic centimeter and temperature of about 3000K, within the range 
well-studied by solar physics). Figure 03] shows the prediction for the CMB power spectrum 
from the standard cosmological model to the power spectrum measured by WMAP5. The 
theoretical curve shown has only six free parameters defining its shape, and matches the 
data astonishingly well. 

The six parameters describing the "vanilla" standard cosmological model are the baryon 
and cold dark matter fractions scaled by /i^, i.e., Qbh"^ and ^ch'^, the dark energy fraction 
JIa; the power law parameters for the power spectrum of initial fluctuations As and Ug, 
and a nuisance parameter corresponding to the optical depth r looking back through our 
universe to the point when the hydrogen atoms were reionized at a redshift of z ~ 11. In 
this simplest possible model that can match the data, neutrinos are assumed to be massless, 
dark energy is assumed to be a cosmological constant, our universe is assumed to be flat, 
and a handful of other assumptions are made as well. The measured shape of the CMB 
power spectrum provides tight constraints on the values of these parameters and also can 
be used to place limits on deviations from the vanilla model assumptions, e.g., massive 
neutrinos or dark energy with w ^ —1. Furthermore, measuring the polarization of the 
CMB provides another complementary source of information. (See lde Oliveira-Costal 12005 
for more details.) 



4.2.4 Other probes 

There are several other important cosmological probes that deserve mention here for com- 
pleteness, although we will not focus on them in depth. 

Perhaps the most important for the discovery of dark energy is the use of type la 
supernovae to trace the expansion of the universe. These stellar explosions are the result of a 
white dwarf star in a binary system accreting mass from its companion. Eventually the white 
dwarf will reach a critical mass, known as the Chandrasekhar limit, where that electron 
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degeneracy pressure can no longer support the stellar material. This triggers an extremely 
luminous explosion that can be seen at cosmological distances. Because the physics of these 
systems dictates that the explosion occurs at a particular mass, the resulting supernovae 
have sufficiently similar properties that they can be used as standard (or standardizable) 
candles: since we know their absolute brightness, measuring their apparent brightness tells 
us the distance. Thus by also measuring their redshift we can do zeroth order cosmology and 
effectively measure the scale factor a (t) as a function of cosmic time. Cosmologists set out 
to do this expecting to measure the deceleration of the expansion rate due to gravity, bu t 
instead found that the expansion was accelerating ( Riess et al. 1998 : Perlmutter et al. 19991 ). 
This discovery provided the first evidence for dark energy and has induced a paradigm shift 
in cosmological thinking over the past decade. 

Another key piece of evidence for the standard cosmological model comes from Big 
Bang nucleosynthesis studies. Our universe's hot dense early phase reached temperatures 
high enough to induce nuclear fusion of light elements, including deuterium, helium, and 
lithium. The fact that our universe is roughly 25% helium-4 by mass, too much to have 
been fused within stars, is by itself evidence that our universe was, in fact, once hot and 
dense. Furthermore, careful measurements of the abundances of deuterium, lithium, and 
helium-3 yield an estimate of the baryonic matt er density that is both independen t of the 
CMB results and in good agreement with them (jTytler et al.lll996l : ISteigma 

Another class of cosmological probes focuses on tracing the present day (or low-redshift) 
matter distribution by means other than the galaxy redshift surveys discussed in §4.2.21 
These include galaxy cluster surveys, gravitational lensing, and Lyman alpha absorption by 
cosmic gas clouds. 

Galaxy cluster surveys aim to detect galaxy clusters - the largest gravitationally bound 
and equilibrated clumps of matter in our universe - via the X-ra.v emission of their hot 
gas ( Haiman et al. 20051 ) their effect on the CMB ( Hu et al. 2007 ). or other means, and 
measure their abundance as a function of mass. This provides a sensitive probe of nonlinear 
structure formation by is complicated by the difficulty of linking the cluster mass to its 
observable properties. 

Gravitational lensing takes advantage of the general relativistic effect that mass bends 
light and aims to map out the mass distribution directly by studying the distortion of the 
light from distant objects. This has the advantage of dodging the issue of galaxy bias that 
plagues the interpretation of galaxy surveys, but is subject to other challenging systematic 



effects ( Benabed and van Waerbekel l2004l : iMandelbaum et all l2005l ^ . 



Measurements of Lyman alpha absorption study the light from distant AGNs that has 
passed through clumps of hydrogen along the way. Hydrogen at different redshifts will ab- 
sorb light at the wavelength of the Lyman alpha spectral line of atomic hydrogen, creating a 
series of absorption lines known as the Lyman alpha forest which provides a one-dimensional 
map of the gas density along the line of sight. This probe of large-scale structure is sen- 
sitive to a somewhat higher redshift than current galaxy surveys, but it is constrained by 
the geometry of requiring an AGN backlight w hich makes it difficult to build up a true 
three-dimensional picture ( Weinberg et al. 20031 ). 

Finally, we mention two additional probes that are likely to be important in future 
cosmological pursuits: 21cm studies of the epoch of reionization and measuring baryon 
oscillations in the galaxy distribution. The epoch of reionization probes the time period 
partway between the time when the CMB was emitted and the time at which we observe the 
most distant astronomical objects. This era is known as the cosmological Dark Ages due 
to the lack of information we have about it. However, we do know that our universe went 
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through a phase of reionization when the first stars and galaxies began to shine, since most 
of the hydrogen is ionized today. There are several experiments currently in the construction 
or planning phases aiming to probe this era by d etecting the 21cm emiss i on line of neu- 
tral hydrogen redshifted to low radio frequencies ( Terzian and Laziol 20061 : Bowman et al 



20071 : iFalcke et al.ll2007l ^. These experiments will open a new window on our universe for 
cosmological studies. 

We conclude our whirlwind tour of other cosmological probes by discussing baryon 
oscillations. The physics behind this probe is as follows: The overall matter distribution is 
determined primarily by dark matter, but the sound wave patterns of the baryon-photon 
plasma observed in the CMB leave an imprint of tiny wiggles in the matter power spectrum 
as well. These wiggles, called baryon oscillations, encode a characteristic length scale: the 
distance a sound wave traveled from the time of the Big Bang until decoupling. This 
length scale can be used as a standard ruler: since we know the physical size, measuring 
its apparent angular size on the sky tells us its distance from us. By measuring this length 
scale as a function of redshift using galaxy surveys, we can map out the expansion history 
of the universe and use it as a probe of zeroth order cosmology, which will hopefully provide 
insights about th e nature of dark energy. Current measurem ents of the baryon oscillations 



at two redshifts (lEisenstein et al.l l2005l : iPercival et al 



20071) have already proved to be a 



quite valuable constraint on cosmological parameters ( Komatsu et al. 20081 ). However, the 
effects of nonlinear matter clustering and galaxy bias discussed in^ 4.1.51 complicate efforts 
to measure this standard ruler to high precision (ISmith et al.ll2007l . l2008l ). Detailed studies 



of galaxy bias such as the one presented in chapter [6] can therefore help us more fully utilize 
this promising probe in the future. 



4.3 Example: combining WMAP and SDSS 

Combining CMB and galaxy survey data is a particularly powerful strategy: the early- 
universe snapshot of the CMB and the present-universe snapshot of the galaxy distribution 
provide complementary information that allows us to precisely constrain the parameters of 
the cosmological model describing the evolution from one state to the other. In this section 
we discuss one exampl e of this: the comb ination of CMB data from the WMAP 3 year 
data rel ease (WMAP3: ISoereel et al.l 120071 ') and the S PSS LRG samp l e from Data Release 



4 (DR4; lAdelman-McCarthy et ahlbood ) analyzed in iTegmark et all (|2006l ). The goals of 



this analysis are twofold: first, measure the matter power spectrum P (k) on scales large 
enough to avoid unfounded assumptions regarding nonlinear clustering and galaxy bias, 
and second, combine this P {k) measurement with the CMB power spectrum measured by 
WMAP3 to constrain the parameters of the standard cosmological model. 



4.3.1 Power spectrum estimation 

The LRG galaxy sample has a number of features that make it ideal for meas uring the 
galax y power spectrum. First of all, its effective volume 14fj {k) (as defined by iTegmark 
19971 ) is larger by a factor of six than the SDSS main galaxy sample and over ten times 
larger than that of the 2dFGRS. The error bars on the power spectrum scale roughly as 

1 /2 

AP (fc) oc Feff {ky ' , so the SDSS LRGs yield the smallest error bars. Additionally, since 
the LRGs are selected to be galaxies of a similar physical type, the effects of luminosity- 
and color-dependent bias explored in chapter [6] should not have a significant effect. 



82 



+ 90 



180° 




180° 



-90° 



Figure 4-7 The angular distribution of the SDSS DR4 luminous red galaxies shown in 
Hammer- Aitoff projection in celestial coordinates, with the s even colors/greys indicating the 
seven angular subsamples analyzed in lXegmark et all constructed with the mangle 

software discussed in chapter [5l 



The angular distribution of the LRGs from DR4 is shown in Fig. 14-71 The regions 
of the sky containing the galaxies illustrate the area that has been covered by the SDSS 
spectroscopic survey at the time of DR4 - the task of managing such angular coverage 
masks is discussed in detail in chapter [5l Here we just note that the updated mangle 
software was used to create the masks for the angular subsamples shown in different colors 
in Fig. 14-71 The method used to calculate t he power spectrum is t he Pseudo Kahrunen- 
Loeve (PKL) eigenmode method described in Tegmark et al. ( 2004bl ) . which, among other 
advantages, produces uncorrelated error bars on the power spectrum. In order to speed up 
the computation time, the PKL method was applied separately on 21 subvolumes (the 7 
angular subsets shown in Fig. 14-71 x 3 radial subsets) and then combined with minimum 
variance weighting. Additionally, a global PKL analysis of the full data set was done to 
specifically target length scales larger than the individual subvolumes. 

The measured power spectra for both the LRGs and the main sample galaxies are 
shown in Fig. 14-81 Theoretical predictions are also shown, both for linear theory (a.k.a. 
first order cosmology; see §4.1.4p and nonlinear theory (a.k.a. higher or der cosmology; see 
The particular flavor of nonlinear niqdelin g used here is based on lEisenstein and Hu 




Cole et al.l tOQ^ ): lEisenstein etaP jioO^). The measured galaxy power spectrum 
is modeled as 

(k) = Pdewiggled (k) fe2 l + Qnlfc^ , (4.55) 



Pa 



1 + lAk 



- friBwiggiBH (k) models th e damping of the baryon oscillations on small scales and is given by 
(jEisenstein et al.ll2007l ) 



i'dcwiggicd (k) = Wik)P (k) + [1-W (k)] 

-fnowigglo 



(4.57) 



wher e P (k) is the linear theory power spectrum calculated using CMBFAST (ISeliak and Zaldarriaea 
19961 ). -Pnowiggie (k) is given by the baryon- free fitting formula in Eisenstein and Hu (|l999l )" 
and W (k) = e~^^/^*'> is a weighted averaging factor combining the two. This retains the 
wiggles on large scales and fades them out beginning around k = k^,. The wiggle suppres- 



sion scale fc* is defined by l/o", where a = o'']^^a}/^ (^s/0.6841)^^^, a± and (T|| are given by 



,1/2 
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equations (12) and (13) in Eisenstein et al. ( 2007I ). and Ag is the power law normalization 



for the primordial spectrum. The parameters b and Qni in equation (I4.56|) parameterize 
bias and nonlinear evolution respectively: b is the usual linear bias factor defined in ^4.1.4^ 
and increasing values of Q^i indicate stronger nonlinear evolution. Ag is fixed by the CMB 
observations ~ 6^ characterizes the difference in normalization between what one predicts 
for the matter power spectrum given the CMB data and what is observed for the galaxy 
power spectrum. 

The theoretical curves for linear theory shown as solid lines in Fig. 14-81 are calculated us- 



ing t he best-fit cosmological parameters from the WMAP3 CMB data alone (jSpergel et al 



20071 ) plus a fit to the data for the normalization 5^, giving b = 1.9 for the LRGs and 
b = 1.1 for the main galaxies relative to the predicted matter power spectrum. The 
nonlinear models shown as dashed lines again use the WMAP3 cosmological parameters, 
plus a fit for b and Q^i, yielding Qni = 4.6 for the main galaxies and Qni = 30 for the 
LRGs. The models agree quite well with the data, which is a good indicator that the 
nonlinear modeling prescription outlined above captures the important effects, but also 
indicate that no nlinear evolution is q uite strong for the LRGs and becomes important 
at larger scales. iTegmark et al.l ( 20061 ) perform several tests to ensure that their cosmo- 
logical results are not sensitive to the details of this largely empirically-based nonlinear 
modeling, but it is important to note that if we want to improve the precision with which 
future galaxy surveys can do cosmology, it will be essential to understand bias and non- 



(Smith et al. 


2007: 


Nishimichi et al. 


2007: 


Tinker et al. 


2007: 


Zhene: and Weinberg 


2007; 


Taruva and Hiramatsu 2008: Crocce and Scoccimarro 20081: Smith et al. 20081). 



4.3.2 Cosmological parameter estimation 



With the galaxy power spectrum in hand, the next goal is to combine it with the WMAP3 
CMB power spectrum to constrain cosmological parameters. For simplicity, only the LRG 
power spectrum on scales larger than k < 0.2/i/Mpc are used here. 

The cosmological model is parameterized in terms of 12 standard parameters, plus the 
two nuisance parameters b and Qni from equation (j4.56p : 



P = {^tot,^A,i^b,i^c,i^u,w,As,r,ns,nt,a,T,b,Qni)- 



(4.58) 



Table 14.11 defines these 14 parameters and another 45 that can be derived from them; 
in essence, {Qtot,^A,'^b,^c,'^u,w) define the cosmic matter budget, {As,ns,a,r,nt) spec- 
ify the seed fluctuations and (r, 6, Qni) are nuisance parameters. The "vanilla" model, 
i.e., the minimal model parametrized by {Q\,u;b,u;c, As,ns,T,b,Qni), setting Wi, = a = 
r = nt = 0, Otot = 1 and w = —1, is the smallest subset of our parameters that 
provides a good fit to the WMAP plus LRG data - this is the same vanilla model dis- 
cussed in §4.2.31 with the addition of b and Qni- Constraints on these parameters from the 
WMAP p lus LRG data were cq i nputed using the Monte Carlo Markov Chain (MCMC) 



approach (Me tropolis et al. 1953 : Has tings^^l97C: Gilks et al. 1996 : Christensen et al. 2001 



Lewis and Bridle 2002 : Slosar and Hobson 2003l : IVerde et al.ll2003l ) as implemented in lTegmark et al.l 
(|2004al ). The best-fit parameter values and their uncertainties are shown in Table UTT]- more 
details about the assumptions that have gone into these parameter estima tions as well as 



tests o f how robust these constraints are to these assumptions can be found in lTegmark et al 



(i2006l ) 
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Table 4.1 Cosmologi c al pa rameters measured from WMAP and SDSS LRG data, repro duced 



from Tegmark et al. ( 20061 ). Reference s cited in the table 
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Figu re 4-8 Measured powe r spectra for the full LRG and main galaxy samples, reproduced 
from Tegmark et al. ( 20061 ). The solid curves correspond to t he liri ear theory cosmological 
model fits to WMAP3 alone from Table 5 of ISpergel et alj (|2007l ). normalized to galaxy 
bias b = 1.9 (top) and b = 1.1 (bottom) relative to t he z = matter power. The dashed 
curves include the nonlinear correction of Cole et al. ( 20051 ) for A = 1.4, with = 30 for 
the LRGs and Qni = 4.6 for the main galaxies; see equation ()4.56p . The onset of nonlinear 
corrections is clearly visible for k > 0.09/i/Mpc (vertical line). 



The parameter constraints provided by any given experiment define an allowed region 
within the A/^-dimensional parameter space (where is the number of parameters) that can 
be quite complicated and may have degeneracies between some parameters ~ for example, 
CMB data constrains ujm = ilm^^ quite tightly but gives only weak constraints on and 
h individually. A common way to illustrate this is to choose two parameters and plot the 
constraints in a two-dimensional plane. Such plots illustrate how c ombining differe r it dat a 
sets breaks various parameter degeneracies, as discussed in depth in iTegmark et"all (|2006l ). 



As an example here, we show just one of these plots in Fig. 14-91 that is particularly 
relevant for particle physics: the constraints on the dark matter fraction uja and the fraction 
fiy of the dark matter th at is composed of n eutrinos. This plot shows that the first year 
WMAP data (WMAPl; ISuergel et all l2003l l only constrains lv^ without saying anything 
about neutrinos, and the WMAP3 data indicates that at most about 20% of the dark 
matter mass could be neutrinos. Adding the SDSS LRG data tightens the constraints still 
further, allowing only about 10% of the dark matter to be neutrinos. Assuming that there 
are no sterile neutrinos, this implies a constraint that the mass of any neutrino must be 
less than ~ 0.3 eV. This is comparable to the sensitivity of the KATRIN particle physics 
experiment current ly being conduc ted to measure the mass of the electron neutrino via 
tritium beta decay ( Lobashev 2003 1). Thus by looking at distant galaxies and the radiation 
that fills our universe, we can do some extremely interesting and relevant particle physics. 
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Figure 4-9 95% constraints in the {oj^, fu) plane, reproduced from Tegmark et al. ( 20061 ). 



The large shaded regions are ruled out by WMAPl (red/dark grey) and WMAP3 (or- 
ange/grey) when neutrino mass is added to the 6 vanilla parameters. The yellow/light grey 
region is ruled out when adding SDSS LRG information. The five curves correspond to M,^, 
the sum of the neutrino masses, equaling 1, 2, 3, 4 and 5 eV respectively - barring sterile 
neutrinos, no neutrino can have a mass exceeding ~ Mjy/3 ~ 0.3 eV (95%). 



4.4 Issues with using galaxies as cosmic tracers 

Measuring the distribution of galaxy clustering has proved to be an extremely powerful tool 
for measuring cosmological parameters. However, there are still issues, both practical and 
theoretical, that need to be addressed in order to optimally extract cosmological information 
from the next generation of galaxy surveys. The following two chapters focus on two of 
these issues: chapter [5] addresses the technical matter of managing angular masks defining 
the survey area on the sky, and chapter [6] addresses the more fundamental question of 
understanding how closely galaxies trace the underlying matter distribution. 

4.4.1 Managing angular masks of galaxy surveys 

As galaxy surveys become larger and more complex, keeping track of the completeness, 
magnitude limit, and other survey parameters as a function of direction on the sky becomes 
an increasingly challenging computational task. For example, typical angular masks of 
SDSS contain about = 300,000 distinct spherical polygons. Managing masks with such 
large numbers of polygons becomes intractably slow, particularly for tasks that run in time 
O {N^) with a naive algorithm, such as finding which polygons overlap each other. In 
chapter [5l we present a "divide-and-conquer" solution to this challenge: we first split the 
angular mask into predefined regions called "pixels," such that each polygon is in only one 
pixel, and then perform further computations, such as checking for overlap, on the polygons 
within each pixel separately. This reduces O (A^) tasks to 0{N), and also reduces the 
important task of determining in which polygon(s) a point on the sky lies from O (A) to 
O (1), resulting in significant computational speedup. 
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Additionally, we present a method to efficiently convert any angular mask to and 
from the popular HEALPix format. This method can be generically applied to con- 
vert to and from any desired spherical pixelization. We have implemented these tech- 
niques in a new version of the MANGLE software package, which is freely available at 
http://space.iiiit.edu/home/tegmark/inaiigle/, along with complete documentation and 
example applications. These new methods should prove quite useful to the astronomical 
community, and since MANGLE is a generic tool for managing angular masks on a sphere, 
it has the potential to benefit terrestrial mapmaking applications as well. 

4.4.2 Relative galaxy bias 

Differences in clustering properties between galaxy subpopulations complicate the cosmo- 
logical interpretation of the galaxy power spectrum, but can also provide insights about the 
physics underlying galaxy formation. Chapter [6] details a study of the nature of this relative 
clustering in which we perform a counts-in-cells analysis of galaxies in SDSS to measure the 
relative bias between pairs of galaxy subsamples of different luminosities and colors. We 
use a generalized test to determine if the relative bias between each pair of subsamples 
is consistent with the simplest deterministic linear bias model, and we also use a maximum 
likelihood technique to further understand the nature of the relative bias between each pair. 

We find that the simple, deterministic model is a good fit for the luminosity-dependent 
bias on scales above ~ 2/i~^Mpc, which is good news for using magnitude-limited surveys 
for cosmology. However, the color-dependent bias shows evidence for stochasticity and/or 
non-linearity which increases in strength toward smaller scales, in agreement with previous 
studies of stochastic bias. Also, confirming hints seen in earlier work, the luminosity- 
dependent bias for red galaxies is significantly different from that of blue galaxies: both 
luminous and dim red galaxies have higher bias than moderately bright red galaxies, whereas 
the biasing of blue galaxies is not strongly luminosity-dependent. These results can be used 
to constrain galaxy formation models and also to quantify how the color and luminosity 
selection of a galaxy survey can impact measurements of the cosmological matter power 
spectrum. 
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Chapter 5 



Methods for rapidly processing 
angular masks of next-generation 
galaxy surveys 



This chapter is adapted from the paper "Methods for rapidly processing an- 
gular masks of next-generation galaxy surveys" by Molly E. C. Swanson, Max 
Tegmark, Andrew J. S. Hamilton, and J. Colin Hill, which has been a ccepted for 



publi cation in the Monthly Notices of the Royal Astronomical Society (jSwanson et al 



5.1 Introduction 



Over the past few decades, galaxy surveys have provided a wealth of information about the 
large-scale structure of our Universe, and the next generation of surveys currently being 
planned promises to provide even more insight. In order to realize the full potential of 
upcoming surveys, it is essential to avoid unnecessary errors and approximations in the way 
they are analyzed. The tremendous volumes of data produced by these new surveys will 
shrink statistical uncertainty to unprecedented levels, and in order to take advantage of 
this we must ensure that the systematic uncertainties can keep pace. The purpose of this 
chapter is to maximize the scientific utility of next-generation surveys by providing methods 
for processing angular masks as rapidly and accurately as possible. 

Angular masks of a galaxy survey are functions of direction on the sky that model the 
survey completeness, magnitude limit, seeing, dust extinction, or other parameters that vary 
across the sky. Th e earliest galaxy red shift surveys - the first Center for Astrophysics red- 
shift survey (C fAl:lHuchra et akHlOsi l and the first Southern Sky Redshift Survey (SSRSl; 



da Costa et al.lll99ll ) - had simple angular masks defined b y boundaries in decl ination and 
Galactic latitude. Th e next generation of surveys - IRAS dStrauss et al.lll992l ^ and PCSz 



(jSaunders et al.1 l200d l - had somewhat more complex masks, with some regions of high 



contamination excluded from the survey. 

The p resent generation o f surv eys - the Two Degree Field Galaxy Reds hift Survey 
fgdFGRS: IColless et aLllioOll . I2OO3I ) and the Sloan Digital Sky Survey (SDSS; lYork et al. 
2003) - consist of photometric surveys that identify galaxies and measure their angular 



positions combined with spectroscopic surveys that measure a redshift for each galaxy to 
determine its distance from us. Angular masks are useful for describing parameters for 



89 



both photometric and spectroscopic surveys - for example, seeing and magnitude hmit are 
key parameters to model in photometric surveys, and the survey completeness - i.e., the 
fraction of photometrically selected target galaxies for which a spectrum has been measured 
- is vital for analyzing spectroscopic surveys. 

The angular masks of SDSS and 2dFGRS consist of circular fields defined by the spec- 
troscopic plates of the redshift survey superimposed on an angular mask of the parent 
photometric survey. 2dFGRS uses the Automatic Plate Measurement (APM) survey (Mad- 
dox et al. 1990a,b, 1996) as its parent photometric survey and covers approximately 1500 
square degrees. The APM angular mask consists of 269 5° x 5° photographic plates which 
are drilled with numerous holes to avoid bright stars, satellite trails, plate defects, and so 
forth. Combining the photographic plates, holes, and spectroscopic fields gives a total of 
3525 polygons that define the spectroscopic angular mask of 2dFGRS. 

The SDSS covers a larger area on the sky - 5740 square degrees in Data Release 5 
(DR%; lAdelman-McCarthv et alll2007l ) - and has a yet more complicated angular mask 
than 2dFGRS. The SDSS photometric survey is done by drift-scanning: each scan across 
the sky covers six long, narrow scanlines, and the gaps between these lines are filled in with 
a second scan slightly offset from the first, producing a "stripe" about 2.5° wide assembled 
from 12 scanlines. In addition to the fairly intricate pattern produced by this scanning 
strategy, there are nearly 250,000 holes masked out of the photometric survey for various 
reasons, plus the circular 3° spectroscopic fields. Combining all of these elements produces 
an angular mask for the spectroscopic survey that contains 340,35 1 polygons. 

To accurately manage the 2dFGRS and SDSS angular masks, [Hamilton and Tegmark 
( 2OO4I ) developed a suite of general-purpose software called mangle, which perfor ms several 
i mport ant procedures on angular masks using computational methods detailed in [Hamilton 
(Il993al lbh. This software has proved to be a valuable resource to the astronomical com- 



munity: it has been used in several analyses of galaxy survey data 1 


Tesmark et al. 


2002. 


2004b|: iMandelbaum et al. 2005: I 


ikasre et al. 2005: Park et al. 2005 


: Conrov et al. 


2005; 


Park et al. 2007: Nishimichi et al.l 


2OO7I: 


Shen et al. 20071: 


Wane et al. 


2007: Tinker et all 



20071 ) . Additionally, it was used extensively in the preparation of the New York University 



Value-Added Galaxy Catalog (NYU-VAGC; Blanton et al. 2005b), which has been used as 
the basis for almost all publications on large-scale structure by the SDSS collaboration. 

However, many of functions in the original version of mangle run in O (A^^) time, which 
becomes quite computationally challenging as the size and complexity of surveys continues 
to increase - computations involving the SDSS mask can take several months of CPU time. 
In this chapter we present new algorithms that can process complicated angular masks such 
as SDSS dramatically faster with no loss of accuracy. Our method is based on splitting an 
angular mask into pixels, reducing the processing time to O (N) by adding an O (N log N) 
preprocessing step. Similar methods based on hierarchical spati al subdivisions have been 
found to be useful in the field of computational geometry (see, e.g., Goodman and O 'Rourkd 
20041 ). but have not previously been applied to angular masks in an astronomy context. 



This w ill be especially useful for next-generation surveys, such as the Dark Energy Sur- 
vey (DES: lDark Energv Survev Collaborationll2005l') and surv e ys done with the Wide- Field 
Multi-Object Spectrograph fWFMQS: lYamamoto et ahl bood: Iciazebrook et a l."2006')._the 



Panoramic Survey Telescope and Rapid Respons e System (P a n-STARRS; 



Kaiser 200 



and t he Large Synoptic Survey Telescope (LSST: ITvsct] I2OO2I : IStubbs et al 
20061 : IlSST Collaborati(^l2006l ). DES, Pan-STARRS, and LSST wih perform photometric 



20041 : iTvson 



surveys and use techniques for estimating redshifts based on the photometric information; 
WFMOS will perform spectroscopic surveys using one of the upcoming photometric sur- 
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number of polygons in mask number of polygons in mask 



Figure 5-1 Speed trials for a series of portions of the SDSS DR5 mask with and without 
pixelization. (a) Time required for pixelization, snapping, balkanization and unification 
of the mask, (b) Time required to identify in which polygon each of the ~400,000 SDSS 
DR5 galaxies lies. Each set of trials is fitted with a power law to show how the processing 
time scales with the number of polygons A^. Also shown on the x-axis are the number of 
polygons in the 2dFGRS mask, the SDSS DR5 mask, and conservative estimates for the 
Pan-STARRS and LSST large-scale structure masks based on scaling up 2dFGRS. 



veys for target selection. The methods we present here are useful for both photometric and 
spectroscopic surveys - keeping track of factors such as seeing and dust extinction could 
prove particularly i mportant for photo metric redshift determinations (jCollister et al.ll2007l : 



Ovaizu et al.ll2008l : iBanerii et al.ll2007l l. 



The proposed large-scale structure survey to be produced by Pan-STARRS will cover 
~30,000 square degrees in 5 wavelength bands, with each field being observed ~50 times 
such that the images can be co-added. A naive scaling up of the 2dFGRS area and number of 
polygons gives an estimate of ~2 x 10^ polygons for the final Pan-STARRS mask. Similarly, 
the LSST large-scale structure survey will cover ~20,000 square degrees in 6 bands, with 
~200 co-added images, suggesting ~8 x 10^ polygons for the LSST mask. The need for 
an improvement in mask processing speed is clearly illustrated in Fig. 15-lt with the old 
algorithms, the projected processing time for the LSST mask would be over 6000 years. 
With our new method, this time is reduced by a factor of ~24,000 to just ten days. 

In addition to developing faster algorithms for processing angular masks, we have also 
integrated the MANGLE utilities w ith HEALPix, a widely used tool for discretizing the 
celestial sphere (jCorski et al.ll2005l ). The methods used by mangle are complementary to 
HEALPix: mangle is best used for functions that are piecewise-constant in distinct regions 
of the sky, such as the completeness of a galaxy survey. In contrast, HEALPix is optimal 
for describing functions that are continuously varying across the sky, such as the cosmic 
microwave background (CMB) or the amount of extinction due to galactic dust. The ability 
to convert rapidly between these two formats allows for easy comparison of these two types 
of data without the unnecessary approximation inherent in discretizing an angular mask. 
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Furthermore, converting a mask into HEALPix format allows users to take advantage of 
pre-existing HEALPix tools for rapidly computing spherical harmonics. 

The spectacular surveys on the horizon are preparing to generate massive, powerful 
datasets that will be made publicly available - this in turn necessitates powerful and intuitive 
general-purpose tools that assist the community to do science with this avalanche of data. 
We provide a such tools with this new generation of the mangle software and describe 
these new tools here. However, this chapter is not a software manual (a manual is provided 
on the MANGLE website) but rather a description of the underlying algori thms. These tools 
have been utilized in recent analyses of SDSS data ( Tegmark et al.l 20061 . discussed in ^4.31 
Swanson et al.l 20081 . discussed in Chapter [6j) ; we are now making them public so others can 



use them as well. 

The outline of this chapter is as follows: in §5.21 we give an overview of the terminology 
we use to describe angular masks and the basic tasks we wish to perform, and in ^5.31 we 
detail our algorithms for accelerating these tasks and quantify their speed. We describe our 
methods for integrating MANGLE with HEALPix in §5.41 and summarize in §5.51 



5.2 Mangle terminology 



The process of defining an angular mask of a galaxy survey in a generic way requires a set 
of standardized terminology. We use the terminology from [Hamilton and Tegmark 
and present a summary of it here. Our formal definition of an angular mask is a union of 
an arbitrary number of weighted angular regions bounded by arbitrary numbers of edges. 
The restrictions on the mask are 



1. that each edge must be part of some circle on the sphere (but not necessarily a great 
circle), and 

2. that the weight within each subregion of the mask must be constant. 



This definition does not cover every theoretical possibility of how a piecewise-constant 
function on a sphere could be defined, but in practice it is sufficiently broad to accommodate 
the design of essentially any galaxy survey. Furthermore, as we discuss in detail in §5.4^ a 
curvilinear angular region (such as a HEALPix pixel) can be well-approximated by segments 
of circles at high resolution. As an example of a typic al angular mask, we show a port ion 
of the SDSS angular mask from Data Release 5 (DR5; lAdelman-McCarthv et al.l 120071 ) in 
Fig. [531 

The fundamental building block of an angular mask is the spherical polygon, which is 
defined as a region bounded by edges that are part of a circle on the sphere. An angular 
mask is thus the union of arbitrarily weighted non-overlapping polygons. For convenience, 
we provide an updated version of a table from [Hamilton and TegmarkI (j2004l ) in Table 15.11 
containing the definitions of key terms used in this chapter. 

The basic procedure used by the mangle software to pr ocess an angular mask consists o f 
the following steps, which are described in greater detail in Hamilton and Tegmark ( 20041 ) : 



1. Snap 

2. Balkanize 

3. Weight 
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Right Ascension Riglit Ascension 



Figure 5-2 A portion of th e SDSS DR5 angular mask (Blanton et al. 2005b; 
Adelman-McCarthv et al. 20071 ). Left: Polygons defining the mask: spectroscopic plates 



and lines delineating different scans and spectroscopic plates are shown in black, and holes 
in the mask are shown in blue/gray. Right: Processed version of the mask, shaded according 
to survey completeness. 



Table 5.1 Definitions of Terms, in Alphabetical Order 
Term Definition 



boundary A set of edges bounding a polygon. 

cap A spherical disk, a region above a circle on the unit sphere, 

circle A line of constant latitude with respect to some arbitrary polar axis on the 

unit sphere. 

edge An edge is part of a circle. A polygon is enclosed by its edges. 

great circle A line of zero latitude with respect to some arbitrary polar axis on the unit 
sphere. A great circle is a circle, but a circle is not necessarily a great circle. 

mask The union of an arbitrary number of weighted polygons. 

pixel A special polygon that specifies some predefined region on the sky as part of 

a scheme for discretizing the unit sphere. Once a mask has been pixelized, 
each polygon in the mask is guaranteed to overlap with exactly one pixel. 

polygon The intersection of an arbitrary number of caps. 

rectangle A special kind of polygon, a rectangular polygon bounded by lines of constant 

longitude and latitude, 
vertex A point of intersection of two circles. A vertex of a polygon is a point where 

two of its edges meet. 

weight The weight assigned to a polygon. The spherical harmonics of a mask are 

the sum of the spherical harmonics of its polygons, each weighted according 
to its weight. A weight of 1 is the usual weight. A weight of signifies an 
empty polygon, a hole. In general the weight may be some arbitrary positive 
or negative real number. 
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Figure 5-3 A cartoon illustrating the process of balkanization (a) with no pixelization, (b) 
with pixelization. 



4. Unify 

The snapping step identifies edges of polygons that are nearly coincident and snaps them 
together so the edges line up exactly. This is necessary because nearly-coincident edges can 
cause significant numerical issues in later computations. In practice, this situation occurs 
when two polygons in a survey are intended to abut perfectly, but are prevented from doing 
so by roundoff errors or numerical imprecision in the mask definition. There are several 
tunable tolerances that control how close two edges can be before they get snapped together 
- these can be adjusted based on how precisely the polygons defining the mask are specified. 

Balkanization is the process of resolving a mask into a set of non-overlapping polygons. 
It checks for overlap between each pair of polygons in the mask, and if the polygons overlap, 
it fragments them into non-overlapping pieces. After this is completed, it identifies polygons 
having disconnected geometry and subdivides them into connected parts. The basic concept 
of balkanization is illustrated in the top panel of Fig. 15-31 The purpose of this procedure 
is to define all of the distinct regions on the sky in which the piecewise-constant function 
we are intending to model might take on a different value. For example, the 2dFGRS and 
SDSS spectroscopic surveys generate masks containing many overlapping circles defining 
each spectroscopic field observed. The SDSS spectrographs can observe 640 objects in each 
field, so if there are more than 640 desired target galaxies in the field, they might not all 
be observed. For example, one field may have spectra for 80 per cent of the targets, and 
a neighboring field may have 90 per cent, but in the region where they overlap all of the 
targets may have been observed. This is how the survey completeness is determined and 
illustrates why balkanization is necessary. 

After the mask is balkanized, weights are assigned to each polygon, representing the 
value of the survey completeness (or any other desired parameter) in that region. The way 
this is done depends on how this information is provided for a given survey. For example, 
the 2dFGRS mask software by Peder Norberg and Shaun Colqj provides a function that 

^ |http : / /magniim . anu . edu . au/~TDFgg/Public/Release/Masks , 
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takes an angular position on the sky and returns the completeness, the magnit ude limit, the 



phot ographic plate number, and the value of the parameter ^ (described in IColless et al 



at that location. This information can be imported into MANGLE by producing a list 
of the midpoints of each of the polygons in the mask, applying the 2dFGRS software to 
calculate the value of the desired function at each midpoint, and assigning these weights to 
the appropriate polygon by using the "weight" routine in mangle. For the SDSS mask, 
the files provided by the NYU-VAGcH (Blanton et al. 2005b) already include weights for 
the survey completeness, so this step is not needed. 

The final step of the processing is unification, which discards polygons with zero weight 
and combines neighboring polygons that have the same weight. While not strictly necessary, 
this procedure clears out unneeded clutter and makes subsequent calculations more efficient. 

After an angular mask has been processed in this fashion, it can be used for function 
evaluation: i.e., given a point on the sky, determine in which one of the non-overlapping 
polygons it lies, and then get the weight of that polygon to obtain the value of the function at 
the input point. It can also be used for creating a random sample of points with the same 
selection function as the survey, calculating Data-Random {DR) and Random-Random 
{RR) angular integrals, and computing the spherical harmonics of the mask. The mangle 
software provides utilities for all of these tasks. 



5.3 Speedup: pixelization 

The tasks of snapping, balkanization, and unification all require comparing pairs of polygons 
in the mask. The brute-force method to accomplish this is simply to compare each polygon 
with every other polygon, which is what the original version of MANGLE did. This naive 
algorithm is O (-/V^) , which is easily sufficient for masks such as the 2dFGRS mask, with 
an AT of a few thousand polygons. However, the SDSS mask has about 100 times as many 
polygons, and points to the need for a cleverer approach. The method we present here is a 
divide-and-conquer approach we dub "pixelization," which processes the mask so that each 
polygon needs to be compared with only a few nearby polygons. 

5.3.1 Pixelization concept 

The underlying concept of pixelization is as follows: before performing any snapping, balka- 
nization, or unification, divide the mask into predefined regions called "pixels" and split 
each polygon along the pixel boundaries such that each polygon is only in one pixel. Then 
for following tasks, polygons need only be compared with other polygons in the same pixel. 

The process of balkanization with and without pixelization is illustrated in Fig. 15-31 for 
three overlapping polygons. The top panel shows the unpixelized version: balkanization 
checks for overlap between each pair of polygons, and then fragments it into seven non- 
overlapping polygons. 

The bottom panel shows the same process but using pixelization as a first step. First, 
each polygon is divided along the pixel boundaries, shown by the dotted lines. At this point, 
each of the four pixels shown has three polygons in it: the intersections between that pixel 
and each of the original three polygons. Then balkanization is performed within each pixel: 
the three polygons in the upper left pixel are split into five non-overlapping polygons, and 
so forth, yielding a final set of 18 non-overlapping polygons. In this illustrative example, 

http: //sdss .physics .nyu. edu/vagc. 
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pixelization increases the complexity of the process, but in general it replaces the O (-/V^) 

algorithm for balkanization with one that is roug hly O {N/Mf^ for N poly gons and 

M pixels, which is roughly O (N) if M ~ A^. For large, complicated masks such as SDSS, 
this speeds up the processing time by a factor of ~1200. 

Once a mask has been pixelized, the important task of determining in which polygon(s) 
a given point lies is sped up greatly as well: one merely has to calculate in which pixel the 
point lies, and then test if the point is in each polygon within that pixel. For a typical 
pixelization scheme, the appropriate pixel number for a given point can be found with a 
simple formula, i.e. an O (1) calculation, so pixelization reduces the O (N) algorithm of 
testing every polygon in the mask to 0{N/M). This means that with pixelization, this 
task does not depend on the total number of polygons in the mask at all if M ~ A^. 

It is important to note that the pixelization procedure makes no approximations to the 
original mask - the pixels are used simply as a tool to determine which polygons are close to 
each other, not as a means of discretizing the mask itself. A mask that has been balkanized 
after pixelization contains all of the same information as it would without pixelization, 
except that it took a tiny fraction of the time to produce. Furthermore, unification can be 
applied across the whole mask rather than within each pixel, which effectively unpixelizes 
the mask if desired. Thus there are essentially no drawbacks to using our pixelization 
procedure. 



5.3.2 Pixelization schemes 



Simple scheme 



The most straightforward means of pixelizing the sky is to use lines of equal azimuth 
and elevation as the pixel boundaries. The azimuth and elevation typically correspond to 
celestial coordinates ~ right ascension (RA) and declination (Dec) - in a survey mask. In 
this scheme, the whole sky is split into quadrants along the equator and prime meridian to 
form the lowest resolution of pixelization. In mangle this is defined as resolution 1. 

To pixelize to higher resolutions, each pixel is split into four child pixels, with the 
boundaries at the midpoints of azimuth and of cos(elevation) within the pixel. This creates 
pixels with equal area. Thus resolution 2 consists of each of the four resolution 1 quadrants 
split into four pixels each, and so forth. This procedure pr oduces a hierarchica l pixel 
structure known in computer science terminology as a quadtree ( de Berg et al. 2000l ) : there 



are 4'' pixels in this scheme at resolution r, and each of these pixels has four child pixels at 
resolution r + 1 and r parent pixels, one at each lower resolution. Resolution is defined 
as being the whole sky. 

In MANGLE, the pixels of the simple scheme are numbered as follows: the whole sky is 
pixel 0, the four quadrants of resolution 1 are pixels 1, 2, 3, and 4, resolution 2 is pixels 5-20, 
and so forth. Thus just one number specifies both the resolution and the pixel location. At 
each resolution, the pixels are numbered in a ring pattern, starting from an elevation of 90° 
and along increasing azimuth for each ring of equal elevation. The simple scheme pixels at 
the four lowest resolutions are shown in Fig. l5-4[ 
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Figure 5-4 The full sky (shown in a Hammer-Aitoff projection in celestial coordinates) 
pixelized with the simple pixelization scheme at the four lowest resolutions. 




Figure 5-5 The full sky (shown in a Hammer-Aitoff projection in celestial coordinates) 
pixelized with the SDSSPix pixelization scheme at the four lowest resolutions. 
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SDSSPix scheme 



Alternatively, the pixelization can be done such that it is more closely aligned with the mask 
of a given survey. In particular, a pixelization scheme called SDSSPi^l has been developed 
for use with the SDSS geometry. Like the simple scheme, SDSSPix is a hierarchical, equal- 
area pixeliz a,tion scheme, an d it is based on the SDSS survey coordinates A and rj. As 
described in IStoughton etHl jiool), SDSS survey coordinates form a spherical coordinate 
system rotated relative to the celestial coordinate system. The poles are located at RA = 
95°, Dec = 0° and RA = 275°, Dec = 0° (J2000), which are strategically located outside 
the SDSS covered area and in the galactic plane, rj is the azimuthal angle, with lines of 
constant 77 being great circles perpendicular to the survey equator, and A is the elevation 
angle, with lines of constant A being small circles parallel to survey equator. A = 0°, 77 = 0° 
is located at RA = 185°, Dec = 32.5° with rj increasing northward. This configuration has 
been chosen such that the stripes produced by the SDSS scanning pattern lie along lines of 
constant rj. 

The SDSSPix base resolution is defined by 36 divisions in the rj direction (equally spaced 
in rj) and 13 in the A direction (equally spaced in cos A), for a total of 468 equal-area pixels. 
These divisions are chosen such that at a special resolution level (called the "superpixel" 
resolution), there is exactly one pixel across each SDSS stripe. Finally, as in the simple 
scheme, higher resolutions are achieved by hierarchically subdividing each pixel at a given 
resolution into four smaller pixels. 

SDSSPix has been included in mangle by incorporating several routines from the 
SDSSPix software package available online.^ The numbering scheme in the mangle imple- 
mentation of SDSSPix differs somewhat from the internal SDSSPix numbering: in mangle, 
the entire sky is pixel 0, as in the simple scheme, but it has 117 child pixels (instead of 4), 
numbered from 1 to 117. These pixels, which comprise resolution 1, are not part of the 
official SDSSPix scheme, but are created by combining sets of 4 pixels from what is defined 
as resolution 2 in mangle. The resolution 2 pixels are the official SDSSPix base resolution 
pixels, and are numbered from 118 to 585. Higher resolutions are constructed through the 
standard SDSSPix hierarchical division of each pixel into 4 child pixels. As in the simple 
scheme, the pixel number identifies both the resolution and the pixel position. The super- 
pixel resolution described above is defined as resolution 5 in mangle and contains a total 
of 7488 pixels. 

Mangle typically uses RA and Dec as its internal azimuth and elevation coordinates, 
so the SDSSPix pixels are constructed as rectangles in — A coordinates and then rotated 
into celestial coordinates. The SDSSPix pixels at the four lowest resolutions are shown in 
Fig. EH 



Other schemes 

The implementation of pixelization in MANGLE is designed to be flexible: it is simple for 
users to add their own scheme as well. The pixelization uses only four basic routines: 

1. get_pixel: Given a pixel number, return a polygon representing that pixel. 

2. which_pixel: Given a point on the sky and a resolution, return the pixel number 
containing the point. 

^httpT7/laliniu.phyast .pitt . edu/~scraiitoii/SDSSPix/| 
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3. get_child_pixels: Given a pixel number, return the numbers of its child pixels. 

4. get_parent_pixels: Given a pixel number, return the numbers of its parent pixels. 

Adding a new pixelization scheme simply requires creating appropriate versions of these 
four routines. Note that it also requires that the pixels can be represented as polygons - 
this is not strictly the case for the HEALPix pixels, as discussed further in ^5.4[ 

5.3.3 Pixelization algorithm 

The purpose of pixelization is to speed up the processing of angular masks, which means 
that the pixelization itself must be done with a clever, speedy algorithm or nothing will be 
gained. The naive algorithm is to search through all of the polygons in the mask for those 
that overlap that pixel. This is O {NM) for N polygons and M pixels, and is not sufficient 
for our purposes. 

Our fast pixelization algorithm is a recursive method that takes advantage of the hier- 
archical nature of the pixelization schemes. The method works as follows: 

1. Start with all the mask polygons that are in pixel i. 

2. Create polygons for each child pixel of pixel i at the next resolution level. 

3. Split the mask polygons in pixel i along the child pixel boundaries such that each 
polygon lies within one child pixel. 

4. Repeat steps 1-3 for the mask polygons in each child pixel until desired stopping point 
is reached. 

Starting this with pixel i = 0, i.e., the whole sky, will pixelize the entire mask in O {N log M) 
time. Thus the pixelization does not add too much overhead time to the overall mask 
processing. 

There are two different methods for choosing the desired stopping point of the pixeliza- 
tion. The simplest method is to stop at a fixed resolution, such that the entire mask is 
pixelized with pixels of the same size. The example mask from Fig. 15-21 is shown in the left 
panel of Fig. 15-61 pixelized to a fixed resolution in the simple scheme. 

Alternatively, the stopping condition can be chosen to be a maximum number of poly- 
gons allowed in each pixel: if there are more than A^max polygons in a given pixel, continue 
the recursion and divide those polygons into pixels at the next resolution level. This results 
in an adaptively pixelized mask, where higher resolutions are automatically used in regions 
of the mask that are more complicated. This method is especially useful for masks with 
varying degrees of complexity in different areas. The example mask from Fig. 15-21 is shown 
again in the right panel of Fig. 15-61 pixelized adaptively with A^ax = 30. 

The implementation of pixelization in mangle allows users to choose either of these 
methods and to select values for the fixed resolution level or for Amax- 

5.3.4 Speed trials 

In order to choose optimal default values for the maximum resolution and Amax as well 
as to demonstrate the dramatic improvements in speed that pixelization provides, we have 
conducted a series of speed trials of the new mangle software. 
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Right Ascension Riglit Ascension 

Figure 5-6 Left: portion of SDSS mask from Fig. 15-21 pixelized to resolution 8 (4^ total pixels 
on the sky) with the simple pixelization scheme. Pixel boundaries are shown in red/light 
gray. Right: the same mask pixelized with the simple pixelization scheme using the adaptive 
resolution method with a maximum of 30 polygons per pixel. 



There are two procedures we are interested in optimizing: firstly, basic processing of a 
mask detailed in ^5.21 involving pixelization, snapping, balkanization, and unification, and 
secondly, the use of the final mask to identify in which polygon a given point lies. We carried 
out trials of these two procedures using both the simple and SDSSPix pixelization schemes 
described in ^5.3.2i For each of these schemes, we tested both the fixed and adaptive 
resolution methods for stopping the pixelization algorithm described in §5.3.31 For the 
fixed resolution method, we measured the time for several different values of the maximum 
resolution, and for the adaptive resolution method, we measured the time as a function of 

The results are shown in Fig. l5-7i (Note that the overhead time for reading and writing 
files and doing general setup has been subtracted - the times shown here are just for the 
primary operations.) From the fixed resolution trials, we see that the optimal resolution 
choice for the SDSS mask is the one that has approximately 10^ total pixels on the sky 
for both the simple and SDSSPix schemes. This corresponds to resolution 9 for the simple 
scheme and resolution 6 for SDSSPix. When using the adaptive method, the choice for 
the maximum number of polygons allowed in each pixel that gives the fastest processing 
is -/Vjnax = 40 for the simple scheme and N^s.^ = 46 for SDSSPix. Overall, the fastest 
choice (by a slight margin) for the SDSS mask is using the SDSSPix scheme with adaptive 
resolution. It is also interesting to note that different mangle processes have different 
optimum values - for example, snapping is fastest when there are fewer polygons in each 
pixel compared to balkanization. 

The impact of our pixelization algorithm is most clearly demonstrated by Fig. 15-li This 
shows the processing time and polygon identification time for a series of selected portions 
of the SDSS mask as a function of the total number of polygons, both with and without 
pixelization - again, overhead time has been subtracted here. Pixelization clearly gives an 
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Figure 5-7 Time required to process the full SDSS DR5 mask with different choices of 
pixelization schemes and methods. The "total" curves are the sum of the pixelization, 
snapping, balkanization, and unification curves. Also shown is the time required to identify 
in which polygon each of the ~400,000 SDSS DR5 galaxies lies (polyid), scaled up by a 
factor of 10. The plots on the left show time vs. the number of pixels at a fixed resolution. 
The plots on the right show time vs. A^maxj the maximum number of polygons in allowed 
in each pixel when pixelizing with the adaptive resolution method. 



101 



improvement in speed that becomes increasingly significant for larger numbers of polygons. 

To quantify this, we fit theoretical models to each series of trials. Without pixelization, 
the processing time for snapping, balkanization, and unification is well-fitted by AN"^ where 
A is a free parameter with a best-fitted value of A = 3.2 x 10~^ CPU seconds. With 
pixelization (and using the adaptive resolution method with the optimal choice for A'^max) 
this is reduced substantially. We fit the times using theoretical models based on the formulas 
in ^5.3.1l assuming M ^ N . We fit the pixelization time with i?A^log4 and the snapping, 
balkanization, and unification time with CN where B and C are free parameters. The 
simple scheme gives {B, C) = (0.0054, 4.7 x lO"'') and the SDSSPix scheme gives {B, C) = 
(0.0045, 4.8 X lO"'') (units are all in CPU seconds) - thus for processing the SDSS mask, 
the SDSSPix scheme is slightly faster than the simple scheme. Our fits give C /B = 0.1, so 
the overall processing time scales like O {N + 0.1A^log4 A^). 

These fitted curves are shown in Fig. 15-11 and allow us to extrapolate estimates for 
processing masks with larger numbers of polygons. For the full SDSS DR5 mask containing 
about 300,000 polygons, pixelization reduces the processing time by a factor of ~1200. The 
improvement for future surveys will be even more dramatic: the mask for the co-added 
LSST survey might contain ^l^f polygons - pixelization reduces the processing time by a 
factor of ~24,000. 

The time for identifying in which polygon each SDSS galaxy lies is significantly sped up 
as well: without pixelization it is fitted by DN where D = 0.012 CPU seconds, although the 
scatter is rather large due to dependence on which polygons are used. With pixelization, it 
is well-fitted by a constant (8.2 CPU seconds for the simple scheme, 7.9 CPU seconds for the 
SDSSPix scheme) - thus the time for polygon identification does not depend on how many 
polygons are in the mask if the number of pixels used is chosen to be roughly proportional 
to the number of polygons. For the SDSS mask, the time to identify the polygons for all of 
the SDSS galaxies is reduced by a factor of nearly 500. 



5.4 Unification with HEALPix and other pixelized tools 



HE A LPix (jCorski et al.l l200,4 ICorski et all Il999l ) is a hierarchical, equal-area, isolatitudi- 
nal pixeliza tion scheme for the sphere motivated by computational challenges in analyzing 
CMB data (jOorski et al.lll999l : iBond et al.lll999l ). Its base resolution consists of 12 equal- 
area pixels; to generate higher resolutions, each pixel at a given resolution is hierarchically 
subdivid ed into four smaller pixels, and it represents an interesting class of spherical pro- 
jections ( Calabretta and Roukema 200?! ). Resolution in HEALPix is defined in terms of 
A'side, the number of divisions along the side of a base-resolution pixel required to reach 
the desired resolution. Because of the hierarchical definition of the higher resolution pixels, 
A'side is always a power of 2, and the total number of pixels at a given resolution is 12N^^^^. 

HEALPix is a very useful scheme in that it allows for fast and accurate astrophysi- 
cal comput ations by means of appropriately disc retizing functions on the sphere to high 
resolution ( Wandelt et al. 19981 : Dore et al. 2001). In particular, HEALPix includes rou 
tines for fast computations of spheri cal harmonics ( Hivon et al. 20021 : Wandelt et al. 2001 



Szapudi et al.ll200ll : iDore et al.ll200lh. It is widely u sed in the analysis of cosmic microwave 
background data from WMAP (ISpergel et al.l 120071') and has recently been used to approx- 
imate galaxy survey masks as well (jPercival et al.l 120071 ). 

Combining HEALPix with mangle is useful because it facilitates comparisons between 
galaxy survey data and the CMB as is done in experiments measuring the integrated Sachs- 
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Sunyaev-Zel'dovich effect (iMyers et al 



Wolfe effect (jPadmanabhan et al.l 12005: iGiannantonio et al.ll2006l: iRassat et al.l 120071 ) . the 



2004 : Reid and Spergel 2006 ). etc., and can also be 



used for generating m asks that block out regions of high dust extinction from galaxy surveys 
( Schlegel et al. 1998 . dust map available in HEALPix format onlin^). 

In general, it can be applied to any task requiring comparison between a piecewise- 
constant function on the sphere to a continuous function sampled on a discretized spherical 
grid. Converting an angular mask into HEALPix format also allows for rapid computations 
of approximate spherical harmonics of the mask using the existing HEALPix tools. 

The implementation of the HEALPix scheme in mangle consists of two components: 

1. A new polygon format, "healpix_weight" , which allows the user to input a list of 
weights corresponding to each HEALPix pixel at a given Nside parameter; 

2. A new utility, "rasterize", which essentially allows the user to pixelize a mask against 
the HEALPix pixels, by means of a different technique than the pixelization method 
described in ^5.31 



Together, these new features allow effective two-way conversion between the HEALPix 
specifications and those of mangle. 



5.4.1 Importing HEALPix maps into mangle 

The structure of the healpix_weight format is quite simple: an input file consists of a list of 
numbers corresponding to the weight of each HE ALPix pixel a t a giv en A'sidc parameter, 
using the nested numbering scheme described in iGorski et al. In addition, the 

definition of the A'^^jde parameter is extended to include 0, which corresponds to a single 
pixel covering the entire sphere, mangle constructs polygons approximately equivalent to 
the HEALPix pixels through the following procedure: 



1. The exact azimuth and elevation of each vertex of a given pixel are calculated using 
the HEALPix utility "pix2vec_nest" ; 

2. The exact azimuth and elevation of each vertex of the four child pixels of the current 
pixel are calculated using the same utility, four of which are the midpoints of the 
edges of the current pixel; 

3. The four vertices of the given pixel are combined with the four midpoints to construct 
the current pixel. Each edge is defined by the circle that pas ses through the two ver- 
tices a nd the midpoint, using the "edges" format described in [Hamilton and Tegmark 
(I200A . 



4. To eliminate spurious antipodal pieces of the polygon defining the pixel, a fifth cap 
is constructed whose axis coordinates are the exact center of the current pixel and 
whose radius is 10~^ radians greater than the distance from the center of the pixel to 
any of the four vertices — this cap thus encloses the entire pixel. 

A similar technique could be applied to incorporate other pixelization schemes not exactly 
described by spherical polygons as well. It is important to note that the pixels constructed in 
this manner are not exactly equivalent to the actual HEALPix pixels - while the hierarchical 



" jhttp : 7/lambda . gsf c . nasa . gov/ product /f oreground/ebv.map . cf iii| 
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and isolatitudinal properties of the HEALPix pixels are preserved, the equal area property 
is not. For example, at A'^^ide = 1, the areas of the approximate pixels differ on average from 
the actual area by about 0.08%. 

However, as the A'side parameter increases, this difference decreases rapidly: at N^i^f, = 
512, the average difference from the actual area is 0.000002%. The boundaries of the both 
the actual HEALPix pixels and our circles approximating them become straight lines in the 
flat-sky approximation, i.e., when the pixel size is much less than 1 radian. Thus for the 
resolutions at which HEALPix is typically used, the difference is totally negligible. 

In applications involving integrations over the sphere, such as calculating spherical har- 
monics, the slight area differences between the pixels can be corrected for in a straight- 
forward manner since the area of each pixel is known (and can be easily extracted using 
mangle's "area" format): simply multiply the value of the function in each pixel by the 
area of that pixel divided by the average pixel area, and then the HEALPix spherical har- 
monics routines will be exact. Furthermore, the paranoid user can get higher precision by 
rasterizing to a higher resolution and then using the "ud_grade" HEALPix utility to obtain 
results for lower resolutions. 

5.4.2 Exporting polygon files as HEALPix maps 

The idea behind the rasterization method used to convert polygon files into HEALPix maps 
is very similar to that of pixelization: to split up the polygons that comprise a mask using 
a given set of pixels, such that afterward each polygon lies in only one pixel. However, 
rasterization is somewhat different, in that afterward the converse statement also holds: 
each pixel contains only one polygon, namely, itself. In particular, rasterization uses an 
arbitrary user-defined spherical pixelization as the pre-determined scheme against which 
to split up the polygons in a given input mask. In general, the user-defined "rasterizer" 
pixels may be from any pixelization of the sphere, but the method was originally developed 
for use with the approximate HEALPix pixels described in the previous section. From a 
practical standpoint, it would be simplest to implement any spherical pixelization that can 
be exactly represented as spherical polygons by following the steps in §5.3.2t however, for 
pixelization schemes that do not possess this property (such as HEALPix), rasterization 
provides an alternate method of implementation. 

The final product of rasterization is a polygon file in which each polygon corresponds to 
one of the rasterizer pixels (either the approximate HEALPix pixels or a user-defined set 
of pixels) and its weight is the area-averaged weight of the input angular mask within that 
pixel. This output can easily be converted into a FITS file read by the HEALPix software 
using a simple script provided with mangle. 

Rasterization consists of the following steps: 

1. Calculate the area of each rasterizer pixel; 

2. Compute the area of the intersection of each input mask polygon with each rasterizer 
(e.g., HEALPix) pixel; 

3. Calculate the area-averaged weight within each rasterizer pixel. 

Note that, as with snapping, balkanization, and unification, this procedure is greatly ac- 
celerated by pixelizing both the rasterizer polygons and the polygons defining the mask to 
the same resolution using one of the pixelization schemes described in §5.3.21 e.g. simple 
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Figure 5-8 Left: Portion of SDSS mask as shown in Fig. l5-2i Right: Portion of SDSS mask 
from Fig. 15-21 as approximated by HEALPix pixels, rasterized to Ngi^ie = 512. 

or SDSSPix. Then step (ii) in the above procedure then involves comparing only polygons 
in the same simple/SDSSPix pixel. Our example mask from Fig. 15-21 (duplicated in the left 
panel of Fig. I5-8P is shown in the right panel of Fig. 15-81 rasterized with HEALPix pixels 
at A'side = 512. The ability to convert between mangle and HEALPix formats allows for 
straightforward comparisons between different types of functions defined on the sphere, e.g. 
angular masks of galaxy surveys and the CMB, as illustrated in Fig. 15-91 

5.5 Summary 

As technologies for surveying the sky continue to improve, the process of managing the 
angular masks of galaxy surveys grows ever more complicated. The primary purpose of this 
chapter has been to present a set of dramatically faster algorithms for completing these 
tasks. 

These algorithms are based on dividing the sky into regions called "pixels" and perform- 
ing key operations only within each pixel rather than across the entire sky. The pixelization 
is based on a hierarchical subdivision of the sky - this produces a quadtree data structure 
that keeps track of which polygons are nearby each other. The preprocessing step of pix- 
elization is O ( A'^ log A'^) for a mask with A^ polygons and ~A pixels, and it reduces the 
mask processing time from O {N'^) to O (A). Furthermore, it reduces the time required to 
locate a point within a polygon from O (A) to O (1). 

This method is exact, i.e., it does not make a discrete approximation to the mask, and 
it takes only a tiny fraction of the computation time. It accelerates the processing of the 
SDSS mask by a factor of about 1200 and reduces the time to locate all of the galaxies 
in the SDSS mask by a factor of nearly 500. It will provide even more dramatic gains for 
future surveys: processing time for the LSST large-scale structure mask could be reduced 
by a factor of about 24,000. 

We have also described a method for converting between masks described by spherical 
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Figure 5-9 The new version of mangle can output angular masks in HEALPix format, 
which ahows for easy comparisons to CMB and other sky map data and allows users to take 
advantage of existing HEALPix tools. Top: The SDSS DR5 completeness mask, raster- 
ized and plotted using HEALPix routin es. Middle: The final 2dFGRS completeness mask 
as determined in Hamilton et al. ( 20081 ). rasterized and plotted using HEALPix routines. 
Bott om: CMB temperat ure difference map measured by WMAP channel 4, with units in 
mK. ISpergel et al.l (|2007l ) All three maps are shown at Ngid^ = 512 in celestial coordinates. 
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polygons and sky maps in the HEALPix format commonly used by CMB and large-scale 
structure experiments. This provides a convenient way to work with both piecewise-constant 
functions on the sky such as the completeness of a galaxy survey and continuously varying 
sky maps such as the CMB temperature. Converting angular masks into HEALPix for- 
mat also allows users to take advantage of existing HEALPix tools for rapidly computing 
spherical harmonics. 

All of the new algorithms and features detailed here have been integrated into the man- 
gle software suite, which is available for free download at http : / / space . mit . edu/home/te gmark/mangle/[ 
This updated software package should prove increasingly useful in the coming years, espe- 
cially as next-generation surveys such as DES, WFMOS, Pan-STARRS, and LSST get 
underway. 
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Chapter 6 



SDSS galaxy clustering: luminosity 
&; color dependence and 
stochasticity 



This chapter is adapted from the paper "SDSS galaxy clustering: luminosity and 
colour dependence and stochasticity" by Molly E. C. Swanson, Max Tegmark, 
Michael Blanton, and Idit Zehavi, which was previously published in the Monthly 



Notic es of the Royal Astronomical Society 385, pp. 1635-1655 ([Swanson et al. 
20081 ) ■ 



6.1 Introduction 



In order to use galaxy surveys to study the large-scale distribution of matter, the relation 
between the galaxies and the underlying matter - known as the galaxy bias - must be 
understood. Developing a detailed understanding of this bias is important for two reasons: 
bias is a key systematic uncertainty in the inference of cosmological parameters from galaxy 
surveys, and it also has implications for galaxy formation theory. 

Since it is difficult to measure the dark matter distribution directly, we can gain insight 
by studying relative bias, i.e., the relation between the spatial distributions of different 
galaxy subpopulations. There is a rich body of literature on this subject tracing back 



many decades (see , e.g., Hubble and Humason 1931 • Davis and Geller 19761 : Hamilton 19881 : 



White et al.lll988l : IPark et al.lll994l : iLovedav et al.lll995l). and been studied extensively in 



recent years as we l l, both theoretically (|Seliakll200ll : Ivan den Bosch et al.ll2003l : ICooravllioOsI : 
Sheth et al.l l200,4 iTinker et al 1 120071 ) and observationally. Such studies have established 
that biasing depends on the type of galaxy under considera tion - for example, early-type, red 



galax i es are more clustered than late-type, blue galaxies (iGuzzo et al.lll997l: iNorberg et al 



20021 : iMadgwick et al.1 l2003l : IConwav et al.1 l200,5l : iLi et al 



2006; 



Croton et al.ll2007ll. and 



17 , f-j — — , I —,,[ s^^; — =^ , — 

luminous galaxies are more clustered than dirn galaxies ( Willmer et al.l 19981: Nor berg et al 



200ll : iTegmark et "al1l2004bl : IZehavi et al1l200,5l : ISeljak et al.ll200,4ISkibba et al.ll2006l ). Since 



different types of galaxies do not exactly trace each other, it is thus impossible for them all 
to be exact tracers of the underlying matter distribution. 

More quantitatively, the luminosity depende nce of bias has be e n measured in the 2 



Degr ee Field Galaxy Redshift Survey (2dFGRS; |Collessetal.''200l') (iNorberg et al. 



2001 



20021 ) and in the Sloan Digital Sky Survey fSDSS: lYork et al.li2000. : .Stoughton et al.ll2002l ) 
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(jTegmark et al.ll2004bl : IZehavi et allboosi : iLi et alJliooi ^ as well as other surveys, and it is 
generally found that luminous galaxies are more strongly biased, with the difference becom- 
ing more pronounce d above L^,, the characteristic luminosity of a galaxy in the Schechter 
luminosity function ( Schechter 19761 ). 

These most recent studies measured the bias from ratios of correlation functions or 
power spectra. The variances of clustering estimators like correlation functions and power 
spectra are well-known to be the sum of two physically separate contributions: Poisson 
shot noise (due to the sampling of the underlying continuous density field with a finite 
number of galaxies) and sample variance (due to the fact that only a finite spatial volume 
is probed). On the large scales most relevant to cosmological parameter studies, sample 
variance dominated the aforementioned 2dFGRS and SDSS measurements, and therefore 
dominated the error bars on the inferred bias. 

This sample variance is easy to understand: if the power spectrum of distant luminous 
galaxies is measured to be different than that of nearby dim galaxies, then part of this mea- 
sured bias could be due to the nearby region happening to be more/less clumpy than the 
distant one. In this chapter, we will eliminate this annoying sample variance by comparing 
ho w different galaxies clus t er in th e same r e gion o f spac e, extending the co unts-in-cells work 
of iTegmark and BroinlevI (Il999ll. 'Blantonl ( 200d'). and IWild et alJ (I2OO5II and the correla- 
tion function w ork of Norberg et al.i ()200ll ). Norberg et al. ( 2002 ) lzehavi et ah ( 20051 ) . and 
Li et all (|200fil ). Here we use the counts-in-cells technique: we divide the survey volume 
into roughly cubical cells and compare the number of galaxies of each type within each cell. 
This yields a local, point-by-point measure of the relative bias rather than a global one as 
in the correlation function method. In other words, by comparing two galaxy density fields 
directly in real space, including the phase information that correlation function and power 
spectrum estimators discard, we are able to provide substantially sharper bias c onstraints . 

This local approach al s o enables us to quant i fy so-called stoc hastic bias ( Pen 19981 : 
Tegmark and PeebleslfliisI : IPekel and Lahavlll999l : lMatsuba^ll999l ). It is well-known that 



the relation between galaxies and dark matter or between two different types of galax- 
ies is not necessarily deterministic - galaxy formation processes that depend on v a riable s 
other than the local matter d ensity give rise to s tocha stic bi as as descri b ed in Pen ( 19981 ). 
Tegmark and Peebles! (| 19981 ). IPekel and Lahavl (jl999l ). and iMatsubl^^ (Il999l ). Evidence 



for stochasti city in the r elativ e bias between e arly- t y pe and late- type galaxies has been 



presented in JWild et al.l tOQ^). IConwav et al.l tOQ^ ). iTegmark and Bromley! (|l999l ). and 

3 



BlantonI (|2000|). Additionally, 



Simon et al 



(|2007l ) finds evidence for stochastic bias between 
galaxies and dar k matter via weak lens i ng. T he time evolution of such stoc hastic bias ha s 
been modeled in Tegmark and Peebles! ( 1998! ) and was recently updated in Simon (|2005! l. 
Stochasticity is even predicted in the relative bias between vir ialized clumps of dark mat' 



ter (halos) and the line arly-evolved dark matter distribution ( Casas-Miranda et al. 2002! : 



Seljak and Warren! 12004 ) . Here we aim to test whether stochasticity is necessary for mod- 



eling the luminosity-dependent or the color-dependent relative bias. 

In this chapter, we study the relative bias as a function of scale using a simple stochastic 
biasing model by comparing pairs of SDSS galaxy subsamples in cells of varying size. Such 
a study is timely for two reasons. First of all, the galaxy power spectrum has recently 
been measured to high prec i sion on large scal e s with the goal of constra ining cosmology 



(jTegmark et al! l2004bl . bood : iBlake et al] boO?! : IPadmanabhan et al.! boO?! ) . As techniques 



continue to improve and survey volumes continue to grow, it is necessary to reduce system- 
atic uncertainties in order to keep pace with shrinking statistical uncertainties. A detailed 
understanding of complications due to the dependence of galaxy bias on scale, luminosity. 
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and color will be essential for m aking precise cosmolqgical inferences with the next gener- 
ation , of galaxy redshift surveys ( Percival et a l. '2004'; lAbazaiian et al 20051 : Percival et al 
20071 : IZheng and Weinberg! l2007l : iMoller et al.l.2007. : .Kristiansen et al. 



20071 ). 



Secondly, a great deal of theoretical progress on models of galaxy formation has been 
made in recent years, and 2dFGRS and SDSS contain a large enough sample of galaxies 
that we can now begin to place robust and det ailed observa tion al constraints on these rn od- 
els. The framework known as the halo model ( Seljak 2000l ) (see Cooray and Sheth 20021 for 
a comprehensive review) provides the tools needed to make comparisons between theory 
and observations. The halo model assumes that all galaxies form in dark matter halos, 
so the galaxy dis tribution can be modeled by first de t ermining the halo distribution - ei- 
ther analyti cally (ICatelan et al.lll998l: [McDonald' '2006'; Smith et ajj 120071 ) or using A^-body 
simulations ( Smith et al. 20031 : Kravtsov et al.. 2004 ; Kravtsov 20o'6l ) - and then populating 
the halos with galaxies. This second step can be done using semi-analytica.1 galaxy for - 
mation models (|Somerville et al.ll200ll : iBerlind et al]l2003l : ICroton et allboOfil : iBaughlboOfil ) 
or with a statistical approach using a model for the halo occupation distribution (HO D) 
(jPeacock and Smithllioool : iBerlind and Weinberg] l2002l: ISefusatti and Scoccimarrd I2OO.5I ) or 
conditional luminosity function (CLF) ( Yang et al. 20031 '" van den Bosch et al. 20031 ) which 
prescribes how galaxies populate halos. 



successfully i n a number of different contexts (IScrantonI 12003 



Gao and White 


2007). it has been applied 


(Scranton 2003 


Collister and Lahav 


2OO5I: 



Tinker et all I2OO6I : ISkibba et all I2OO6I ). The correlation between a galaxy's e nvironment 



(i.e., the local density of surrounding galaxies) and its color and luminosity (IHogg et al 



2004 : Blanton et al. 2005a) has been interpreted in t he cont ext of the halo model (iBerliiid et al 



2OO.5I: iBlanton et al.l lioOfl : lAbbas and ShethI lioO^ ). and Ivan den Bosch et all (I2OO3I ) and 



Cooravl (|2005l )" make predictions for the bias as a func t ion o f galaxy type and luminosity 



using the C LF formalism. A d dition ally, Zehavi et al. ( 20051 ). Magliocchetti and Porciani 



(2003 ). and lAbbas and ShethI (|200,5l ) use correlation function methods to study the lumi- 
nosity and color dependence of galaxy clustering, and interpret the results using the Halo 
Occupation Distribution (HOD) framework. The analysis presented here is complementary 
to this body of work in that the counts-in-cells method is sensitive to larger scales, uses a 
different set of assumptions, and compares the two density fields directly in each cell rather 
than comparing ratios of their second moments. The halo model provides a natural frame- 
work in which to interpret the luminosity and color dependence of galaxy biasing statistics 
we measure here. 

The rest of this chapter is organized as follows: ^6.21 describes our galaxy data, and 
§6.3.11 and §6.3.21 describe the construction of our galaxy samples and the partition of 
the survey volume into cells. In ^6.3.31 we outline our relative bias framework, and in 
§6.3.41 and §6.3.5l we describe our two main analysis methods. We present our results in §6.41 
and conclude with a qualitative interpretation of our results in the halo model context in 

MM 



6.2 SDSS Galaxy Data 



The SDSS dYork et al.ll200nl : IStoughton et al.i 
19981 ) on a dedicated telescope dCun n et al 

bandpasses denoted u, g, r, i and z (iFukugita et al. 19961 ). After astrometric calibration 



2002 ) uses a mosaic CCD camera ( Gunn et al 



20061) to ima ge the sky in five photometric 
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(IPier et al.l l2003ll . photometric data reduction (ILupton et al.ll200ll'), and photom etric cah- 
bration (|Hogg et alJlioOll : ISmith et al]l2nn2l : llvezic et alJl2nn4l : IXucker et al.lbnnfil ) , galaxies 
are selected for spectroscopic observations. To a good approximation, the main galaxy 
sample consists of all galaxies with r-band appareii t Petrosian magnitude r < 17.77 after 
correction for reddening as per Schlegel et al. ( 19981 ): there are about 90 such galaxies per 
square degree, with a median redshift of 0.1 and a tail out to z ^ 0.25. Gala xy spectra are 



also measured for the Luminous Red Galaxy sample faisenstein et ahl B. which is not 
used in this chapter. These targets are assigned to spectroscopic plates of diameter 2.98° by 
an adaptiv e tiling algorithm (B lanton et al. 2003b) and observed with a pair of CCD spec- 
trographs ( Uomoto et al. 20041 ) . after which the spectroscopic data reduction and redshift 
determination are performed by automated pipelines. The rms galaxy redshift errors are of 
order 30kms~^ for main galaxies, hence negligible for the purposes of the present chapter. 

Our analysis is based on 380,614 ma in galaxies (the 'safeO' cut) from the 444,189 galaxies 
in the 5th SDSS data release ('DR5') ( Adelman-McCarthv et aP 2007 ). processed via the 
SDSS data repository at New York University (Blanton et al. 2005b). The details of how 



these s amples were processed and ra odeled are given in Appendix A of iTegmark et al 



(|2004bl ;i and in lEisenstein et al.l (|2005l ). The bottom line is that each sample is completely 
specified by three entities: 

1. The galaxy positions (RA, Dec, and comoving redshift space distance r for each 
galaxy) ; 

2. The radial selection function n(r), which gives the expected (not observed) number 
density of galaxies 3jS Si function of distance ; 

3. The angular selection function n(r), which gives the compl eteness as a function of 
direc tion in the sky, specified in a set of spherical polygons (jHamilton and Tegmark 
20041 '). 



The three-dimensional selection functions of our samples are separable, i.e., simply the 
product n (r) = n (r) n (r) of an angular and a radial part; here r = |r| and r = r/r are the 
radial comoving distance and the unit vector corresponding to the position r. The volume- 
limited samples used in this chapter are constructed so that their radial selection function 
n (r) is constant over a range of r and zero elsewhere. The effective sky area covered is 
17 = J n{r) dQ ^ 5183 square degrees, and the typical completeness n (f) exceeds 90 per 
cent. The conversion from redshift z to comoving distance was made for a flat ACDM 
cosmological model with = 0.25. Additionally, we make a first-order correction for 
redshift space distortio ns by applying the finger-of-god compression algorithm described in 
Tegmark et"all (j2004bl ) with a threshold density of 6c = 200. 



6.3 Analysis Methods 

6.3.1 Overlapping Volume-Limited Samples 

The basic technique used in this chapter is pairwise comparison of the thre e-dimensional 
density fields of galaxy samples with different colors and luminosities. As in lZehavi et al 



diooi), we focus on these two properties (as opposed to morphological type, spectral type, 
or surface brightness) for two reasons: they are straightforward to measure from SDSS data, 
and recent work (Blanton et al. 2005a) has found that luminosity and color is the pair of 
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Comoving distance r (h Mpc) 



Figure 6-1 Histogram of the comoving number density (after finger-of-god compression) of 
the volume-limited samples L1-L7. The cuts used to define these samples are shown in 
Table [UTTl Note that the radial selection function n (r) is uniform over the allowed range for 
each sample. Arrows indicate volumes V1-V6 where neighboring volume-limited samples 
overlap. 



properties that is most predictive of the local overdensity. Since color and spectral type 
are strongly correlated, our st udy of the color dependence of bias probes similar physics 
as studies using spect r al type (Tegrnark and Bromlev 1999 : Blanton 200d : Norberg et al 
200i; IWild et ahlEini : IConwav et al.1 bnO.'^l ) . 



Our base sample of SDSS galaxies ('safeO') has an r-band app a rent rn agnitude range 
of 14.5 < r < 17.5. Following the method used in iTeemark etaP (|2004bl l. we created a 
series of volume-limited samples containing galaxies in different luminosity ranges. These 
samples are defined by selecting a range of absolute magnitude Miuminous < Mo.i^ < M^im 



and defining a redshift range such that the near limit has Mo.i, 



M, 



luminous ) ^ 



14.5 



and the far limit has Mo.: 



Mdi, 



17.5. Thus by discarding all galaxies outside the 



redshift range, we are left with a sample with a uniform radial selection function n (r) that 
contains all of the galaxies in the given absolute magnitude range in the volume defined by 
the redshift limits. Here Mo.ij. is defined as the absolute magnitude in the r-band shifted 
to a redshift of z = 0.1 (Blanton et al. 2003a). 

Our volume-limited samples are labeled LI through L7, with LI being the dimmest and 
L7 being the most luminous. Figure [6^ shows a histogram of the comoving galaxy density 
n (r) for L1-L7. The cuts used to make these samples are shown in Table [6Tl 

Each sample overlaps spatially only with the samples in neighboring luminosity bins - 
since the apparent magnitude range spans three magnitudes and the absolute magnitude 
ranges for each bin span one magnitude, the far redshift limit of a given luminosity bin is 
approximately equal to the near redshift limit of the bin two notches more luminous. (It is 
not precisely equal due to evolution and K-corrections.) 

The regions where neighboring volume-limited samples overlap provide a clean way 
to select data for studying the luminosity-dependent bias. By using only the galaxies 
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Table 6.1 Summary of cuts used to create luminosity-binned volume-limited samples. 



Luminosity-binned 



volume- limited 
samples 


Absolute ma 


2;nitude 


Redshift 


Comoving number 
density n yn^Mpc ) 


LI 


-17 < Mo.ir 


< 


-16 


0.007 <z< 0.016 


(1.63 ±0.05) X 10"^ 


L2 


-18 < Mo.v 


< 


-17 


0.011 < z < 0.026 


(1.50 ±0.03) X 10-2 


L3 


-19 < Mo.i., 


< 


-18 


0.017 < z < 0.041 


(1.23 ±0.01) X 10-2 


L4 


-20 < Mo.v 


< 


-19 


0.027 < z < 0.064 


(8.86 ± 0.05) X 10-3 


L5 


-21 < Mo.v 


< 


-20 


0.042 < z < 0.100 


(5.02 ± 0.02) X 10-3 


L6 


-22 < Mo.ir 


< 


-21 


0.065 < z < 0.152 


(1.089 ±0.005) X 10-3 


L7 


-23 < Mo.i, 


< 


-22 


0.101 < z < 0.226 


(4.60 ± 0.06) X 10-^ 



Table 6.2 Overlapping volumes in which neighboring luminosity bins are compared. 



Pairwise comparison 
(overlapping) volumes 


Overlapping 
bins 


Redshift 


VI 


LI & L2 


0.011 <z< 0.016 


V2 


L2 & L3 


0.017 < z < 0.026 


V3 


L3 k L4 


0.027 < z< 0.041 


V4 


L4 & L5 


0.042 < z < 0.064 


V5 


L5 & L6 


0.065 <z< 0.100 


V6 


L6 & L7 


0.101 < z< 0.152 



in the overlapping region from each of the two neighboring luminosity bins, we obtain 
two sets of objects (one from the dimmer bin and one from more luminous bin) whose 
selection is volume-limited and redshift-independent. Furthermore, since they occupy the 
same volume, they are correlated with the same underlying matter distribution, which 
eliminates uncertainty due to sampl e variance and ren ioves potential systematic effects due 
to sampling different volume sizes (j Jovce et al.l l2005l ) . We label the overlapping volume 
regions VI through V6, where VI is defined as the overlap between LI and L2, and so 
forth. The redshift ranges for V1-V6 are shown Table [6^21 

To study the color dependence of the bias, we further divide each sample into red galaxies 
and blue galaxies. Figure [6^ shows the galaxy distribution of our volume-limited samples 
on a color-magnitude diagram. The sharp boundaries between the different horizontal 
slices are due to the differences in density and total volume sampled in each luminosity 
bin. This diagram illustrates the well-known color bimodality, with the re dder galaxies 



fallirig predominantly iii a region coinmonl y known as the E-SO ridgeline (jBalogh et al 



2004ljMateus et ahlljoO^ : iBaldrv et al.ll2006l ^. To separate the E-SO ri dgeline from the res t 



of the population, we use the same magnitude-dependent color cut as in IZehavi et all (|200,5l ): 
we define galaxies with ^-^ {g — r) < 0.9 — 0.03 (Mo.i^ ± 23) to be blue and galaxies on the 
other side of this line to be red. 

In each volume V1-V6, we make four separate pairwise comparisons: luminous galaxies 
vs. dim galaxies, red galaxies vs. blue galaxies, luminous red galaxies vs. dim red galaxies, 
and luminous blue galaxies vs. dim blue galaxies. The luminous vs. dim comparisons 
measures the relative bias between galaxies in neighboring luminosity bins, and from this 
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0.2 0.4 0.6 

0.1 



(g-r) 



Figure 6-2 Color-magnitude diagram showing the number density distribution of the galaxies 
in the volume-limited samples. The shading scale has a square-root stretch, with darker 
areas indicating regions of higher density. The line shows the color cut of (g — r) = 
0.9 — 0.03 (Mo.ij, -|- 23). We refer to galaxies falling to the left of this line as blue and ones 
falling to the right of the line as red. 



we can extract the luminosity dependence of the bias for all galaxies combined and for red 
and blue galaxies separately. The red vs. blue comparison measures the color-dependent 
bias. This set of four different types of pairwise comparisons is illustrated in Fig. 16-31 for 
V4, and the number of galaxies in each sample being compared is shown in Table [ 



6.3.2 Counts-in-Cells Methodology 

To compare the different pairs of galaxy samples, we perform a counts-in-cells analysis: we 
divide each comparison volume into roughly cubical cells and use the number of galaxies 
of each type in each cell as the primary input to our statistical analysis. This method 
is complementary to studies based on the correlation function since it involves point-by- 
point comparison of the two density fields and thus provides a more direct test of the local 
deterministic linear bias hypothesis. We probe scale dependence by varying the size of the 
cells. 

To create our cells, we first divide the sky into two-dimensional 'pixels' at four dif- 
ferent angular resolutions using the SDSSPix pixelization schem^ as implemented by the 
updated version of the angular mask processing so ftware mangle discussed in Chapter [5] 
(jHamilton and Tegmarkll2004l : ISwanson et al.ll2007l ). The angular selection function n (f) is 
averaged over each pixel to obtain the completeness. To reduce the effects of pixels on the 
edge of the survey area or in regions affected by internal holes in the survey, we apply a cut 
on pixel completeness: we only use pixels with a completeness higher than 80 per cent (50 
per cent for the lowest angular resolution). Figure [6^ shows the pixelized SDSS angular 
mask at our four different resolutions, including only the pixels that pass our completeness 



^See http://lahmu.phyast.pitt.edu/ scranton/SDSSPix/ 
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Figure 6-3 Galaxy distributions (after finger-of-god compression) plotted in comoving spa- 
tial coordinates for a radial slice of the volume- limited samples L4 (smaller dots, radial 
boundaries denoted by dashed lines) and L5 (larger dots, radial boundaries denoted by 
solid lines), "which overlap in volume V4. Four different types of pairwise comparisons are 
illustrated: (a) luminous galaxies (L5) vs. dim galaxies (L4), (b) red galaxies vs. blue 
galaxies (both in V4), (c) luminous red galaxies (L5) vs. dim red galaxies (L4), and (d) 
luminous blue galaxies (L5) vs. dim blue galaxies (L4). The shaded regions denote the vol- 
ume in which the t'wo sets of galaxies are compared. A simple visual inspection shows that 
the different samples of galaxies being compared generally appear to cluster in the same 
physical locations - one key question we aim to answer here is if these observed correlations 
can be described with a simple linear bias model. 
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Table 6.3 Number of galaxies in each sample being compared. 





All split by luminosity 


All split by color 




Luminous 


Dim 


Red 


Blue 


V i 


427 


651 


125 


953 


V z 


2102 


2806 


1117 


3791 


Vo 


6124 


8273 


5147 


9250 


V4 


12122 


23534 


17144 


18512 


Vo 


11202 


53410 


37472 


27140 


Vo 


1784 


38920 


27138 


13566 




Red split by luminosity 


Blue split by luminosity 




Red luminous 


Red dim 


Blue luminous 


Blue dim 


VI 


72 


53 


355 


598 


V2 


620 


497 


1482 


2309 


V3 


2797 


2350 


3327 


5923 


V4 


6848 


10296 


5274 


13238 


V5 


7514 


29958 


3688 


23452 


V6 


1451 


25687 


333 


13233 



cut. The different angular resolutions have 15, 33, 157, and 901 of these angular pixels 
respectively. At the lowest resolution, each pixel covers 353 square degrees, and the angu- 
lar area of the pixels decreases by a factor of 1/4 at each resolution level, yielding pixels 
covering 88, 22, and 5 square degrees at the three higher resolutions. 

To produce three-dimensional cells from our pixels, we divide each comparison volume 
into radial shells of equal volume. We choose the number of radial subdivisions at each 
angular resolution in each comparison volume such that our cells are approximately cubical, 
i.e., the radial extent of a cell is approximately equal to its transverse (angular) extent. This 
procedure makes cells that are not quite perfect cubes - there is some slight variation in the 
cell shapes, with cells on the near edge of the volume slightly elongated radially and cells 
on the far edge slightly flattened. We state all of our results as a function of cell size L, 
defined as the cube root of the cell volume. At the lowest resolution, there is just 1 radial 
shell for each volume; at the next resolution, we have 3 radial shells for volumes V4 and V5 
and 2 radial shells for the other volumes. There are 5 radial shells at the second highest 
resolution, and 10 at the highest. 

Since each comparison volume is at a different distance from us, the angular geometry 
gives us cells of different physical size in each of the volumes. At the lowest resolution, where 
there is only one shell in each volume, the cell size is 14/i~^Mpc in VI and 134/i~^Mpc 
in V6. At the highest resolution, the cell size is 1.7/i-iMpc in VI and 16/i-^Mpc in V6. 
Figure 16-51 shows the cells in each volume V1-V6 that are closest to a size of ~ 20/i"^Mpc, 
the range in which the length scales probed by the different volumes overlap. (These are 
the cells used to produce the results shown in Fig. 16-81 ) 

6.3.3 Relative Bias Framework 

Our task is to quantify the relationship between two fractional overdensity fields 5i (x) = 
pi (x) /pi — 1 and 62 (x) = p2 (x) /p2 — I representing two different types of objects. This 
framework is commonly used with types (1,2) representing (dark matter, galaxies), or as in 
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Figure 6-4 The SDSS DR5 angular mask pixelized at the four different resolutions used to 
partition the survey into cells, shown in Hammer- Aitoff projection in equatorial coordinates. 
Shading indicates completeness level: per cent is black, 100 per cent is white. 
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Blantonl ( 200Cll ) . Wild et al. ( 20051 ). and Conway et al. ( 20051 ). (early- type galaxies, late- type 
galaxies). Here we use it to represent (more luminous galaxies, dimmer galaxies) or (red 
galaxies, blue galaxies) to compare the samples described in ^6.3.11 Galaxies are of course 
discrete objects, and as customary, we use the continuous field pa (x) (where a=l or 2) to 
formally refer to the intensity of the Poisson point process involved in distributing the type 
a galaxies, as described in more detail below. 

The simplest (and frequently assumed) relationship between 5i and 82 is linear deter- 
ministic bias: 

62 (x) = h^Ji (x) , (6.1) 

where 5iin is a constant parameter. This model cannot hold in all cases - note that it can 
give negative densities if bw^ > 1 - but is typically a reasonable approximation on cosmo- 
logically large length scales where the density fluctuations (5i <C 1, as is the case for the 
measure ments of the large scale power spectrum recently used to const rain cosmological pa- 
rameters (jTegmark et al.ll2004bl : ICole et al.ll2005l : iTegmark et ahlliooi ). This simple model 
is likely to break down on length scales small enough to be affected by galaxy formation 
physics, and such a breakdown could potentially impact the cosmological interpretation of 
medium- and small-scale power spectrum results from galaxy surveys. 



m 



More complicated mod els allow for non-linearity and stochasticity, as described in detail 
Dekel and Lahavl (|l999l '): 



52(x) =6[(5i (x)]5i (x) + e(x). 



(6.2) 



where the bias h is now a (typically slightly non-linear) function of 5i. The stochasticity is 
represented by a random field e - allowing for stochasticity removes the restriction implied 
by deterministic models that the peaks of 5i and 82 must coincide spatially. Stochasticity 
is basically the scatter in the relationship between the two density fields due to physical 
variables besides the local matter densi ty. Non-loc a l gala xy formation processes can also 
give rise to stochasticity, as discussed in iMatsubaral 



The theoretical model we use to describe the galaxy distribution is based on the idea 
that our universe is one realization drawn from an ensemble of all possible universes that 
have a given set of cosmological parameters. Any measurement we make in our universe can 
be described as a random variable over this ensemble of universes. The value we observe is 
a realization of this random variable. Since we only have one universe to observe, we cannot 
make repeated observations to draw different values from the distribution. However, we can 
make theoretical predictions for the probability distribution for a random variable, and then 
perform statistical tests to see if it is plausible that we would observe the particular value 
we do if our model is correct. This is the basic idea behind the statistical tests detailed in 
gSAland g031 

Our statistical model is composed of two separate random processes. The first compo- 
nent is the overdensity field 6 (x) , which is a realization of a zero- mean random field with a 
power spectrum predicted by a set of cosmological parameters, as described in §4.1.41 Here, 
as in ^4.1.4^ we use regular angle brackets () to denote expectation values over the probabil- 
ity distribution for this random field, e.g., a zero-mean random field 5 (x) has {5 (x)) = 0. 
The second component is a random point process which describes the locations of galaxies. 
This is typically taken to be a Poisson point process, which means that the probability of 
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finding m galaxies in a bounded region A is given by 



P (m galaxies in A) = lAl^^g^^^^) , (6.3) 

ml 

where 

A (^) = / A (x) d^x. (6.4) 
J A 

The function A (x) is the intensity function of the point process. For describing the galaxy 
distribution, we use A (x) = n (x) (1 + 5(x)), where n (x) is the selection function of the 
survey as defined in ^6.21 - in other words, the expected galaxy density if there were no 
clustering. Expectation values over this Poisson point process will be denoted by bold 
angle brackets (). 

Our basic observational measurement is N_a , the number of galaxies of type a in cell 
i observed in our SDSS data sample. We will denote quantities calculated using our SDSS 
counts-in-cells data - i.e., realizations of random variables in our universe - with underlines. 
Non-underlined quantities indicate the random variables themselves, for which we compute 
theoretical expectation values based on the statistical model described above. For example, 
taking the expectation value over the Poisson point process for type a galaxies in cell i gives 

(N!i^)=N(^)(l + 6^) (6.5) 



where Sa^ is the volume average of 6a (x) over cell i and N^^ is the integral of the selection 
function n (x) over cell i. Note that aside from slight variations in the angular selection 

- (i) 

function, Na is the same for all cells i - it is simply the expected number of type a 
galaxies in cell i in the absence of any clustering. Taking the expectation over the random 
field realization as well gives 

^(iV«>)=iVl^). (6.6) 
We estimate the overdensity of galaxies of type a in cell i by 



The n-dimensional vectors 



w ^ . (6.7) 

'a 



(6. 



(n) 



—a 



contain the counts-in-cells data to which we apply the statistical analyses in 
g633]and gHXH 

is the realization of the random variable ga for our particular universe that we 
measure from the SDSS data. The covariance matrix of g^, which will be the same for any 
universe drawn from the ensemble defined by our cosmological model, is given by 

where (5^/3 is a Kronecker delta and is known as the shot noise covariance matrix. The 
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presence of 5ap implies that we are assuming that type 1 galaxies and type 2 galaxies 
are distributed via independent point processes, so their shot noise is uncorrelated - the 
correlations between type 1 and type 2 galaxies are encoded in the relationship between 
between 5i (x) and 82 (x) as in equations (16. ip and (j6.2p . Although one might expect the 
fact that the counts of type 1 and type 2 galaxies in a cell is constrained to be equal to the 
total number of galaxies in the cell could induce correlations in the shot noise, we do not 
explicitly use the combined total count in our analyses - uncorrelated shot noise is thus a 
reasonable assumption. 

For Poisson point processes, we have 

(6.10) 

For comparing pairs of different types of galaxies, we construct the 2n-dimensional data 
vector 




(6.11) 



which, under our theoretical statistical model described above, has covariance matrix 

C^«gg^» = S + N, (6.12) 

with 

and the elements of the matrix S given by 

^l, = (^^^)- (6-14) 

N is a diagonal matrix, which indicates that the shot noise is uncorrelated between different 
cells and also between type 1 and type 2 galaxies within the same cell. 

Regarding the matrix S, other counts-in-cells analyses often assume that the cosmolog- 
ical correlations between different cells can be ignored, i.e., {8'£' 5^p^ = unless i = j. 

Here we account for cosmological correlations by computing the elements of S using the 
best-fitting ACDM matter power spectrum as we will now describe in detail. The power 

spectrum P^^ (k) is defined as (^6a (k) 6^ (k')^) = (27r)^ 6^ (k - k') P^p (k), where 5„ (k) = 
f g-«k xj^ jg ^YiQ Fourier transform of the overdensity field. Pu (k) and P22 (k) are 

the power spectra of type 1 and 2 galaxies respectively, and P12 (k) is the cross spectrum 
between type 1 and 2 galaxies. We assume isotropy and homogeneity, so that Pap (k) is a 
function only of A; = |k|, and rewrite the galaxy power spectra in terms of the matter power 
spectrum P (k): 



V 312 i>22 / 



Puik) = hi{kfP{k) 

Pi2{k) = bi{k)b^{k)n2{k)P{k) 

P22{k) = b2{kfP{k), (6.15) 

which defines the functions 61 (k), 62 (k), and ri2 (k). 

To calculate (^k^Sg^/ exactly, we need to convolve 5a (x) with a filter representing cell 
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i and (5/3 (x) with a filter representing cell j. This is complicated since our cells, while all 
roughly cubical, have slightly different shapes. We therefore use an approximation of a 
spherical top hat smoothing filter with radius R: w{r,R) = 3/{4:TtR^)@{R — r) with the 
Fourier transform given by 



w{k,R) 



{kRf 



[sin(A;i?) - kRcos (kR)] . 



(6.16) 



R is chosen so that the effective scale corresponds to cube s with side len gth L: R = y^5/12L, 
where L is the cell size defined in g6321 (See p. 500 in IPeacockl [l999l for derivation of the 
y5/12 factor.) This gives 



where rij is the distance between the centers of cells i and j. The kernel of this integrand 
- meaning everything besides Pai3 {k) here - typically peaks at /c ~ 1/-R and is only non- 
negligible in a range of Alogj^gfc ~ 1. Assuming that the functions bi{k), b2{k), and 
^12 (k) vary slowly with k over this range, they can be approximated by their values at 
kpeak = 1/R = Y^12/5/L and pulled outside the integral, allowing us to write 



S = (L) 



Sm 

brci (L) rrc\ (L) Sm 



brel (L) Trel (L) S 



M 



brcl{L) S 



M 



(6.18) 



where af (L) = bi (/Cpeak)^, ^'rel (L) = 62 (^peak) /bl (kpenk), rrel (L) 

the correlation matrix for the underlying matter density: 



ri2 (kpesk), and Sm is 



1 

2^ 



sin (kri 
kvi-i 



'-P{k) \w {k,R)\^k'^ dk 



(6.19) 



non-linear transformation of Smith et al. 



For th e matter power spectrum P{k), we use the fi tting formula from Novosvadlvj et al.l 
( I999I ) with the best-fitting ' vanilla' parameters from Tegmark et al. ( 2004a ) and apply the 



In principle, using P (k) as determined 
previously from SDSS data could introduce some model dependence in our results. However, 
our tests in §6. A. II using an uncorrelated (i.e., diagonal) signal matrix indicate that the 
particular choice of the matter power spectrum has no significant effect on the results, so 
this is not a problem in practice. 

Our primary parameters are the relative bias factor 6rei {L) , the relative cross-correlation 
coefficient rrei(-C'), and the overall normalization af{L). The only assumptions we have 
made in defining these parameters are homogeneity, isotropy, and that 61 (k), 62 (k), and 
ri2 (k) vary slowly in k. These parameters are closely related to those in the biasing models 
specified in equations (j6.ip and (j6.2p : If linear deterministic biasing holds, then b^ei = ^Un 
and Trei = 1, and the addition of either non-linearity or stochasticity will give r^^i < 1. As 
discussed in Matsubara ( 19991 ). stochasticity is expected to vanish in Fourier space (i.e., 
^12 {k) = 1) on large scales where the density fluctuations are small, but scale dependence 
of bl (k) and 62 (k) can still give rise to stochasticity in real space. We will measure the 
parameters ?Jrei (L) and r^-^i (L) as a function of scale, thus testing whether the bias is scale 
dependent and determining the range of scales on which linear biasing holds. 
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6.3.4 The Null-buster Test 



Can the relative bias between dim and luminous galaxies or between red and blue galaxies 
be explained by simple linear determ inistic biasing? To address th is question, we use the 
so-called null-buster test described in Tegmark and Bromley ( 19991 ). For a pair of different 



types of galaxies, we calculate an n-dimensional difference map vector 

gA = &-/gi (6.20) 

for a range of values of /. (Note that contains a model parameter / even though it 
is constructed from observed quantities.) If equation (j6.ip holds and / = ftun, then the 
density fluctuations cancel and will contain only shot noise, with a covariance matrix 
((gAgA)) = I^A = N2 + /^Ni - this is our null hypothesis. Alternatively, if equation (|6.1|) 
does not hold, the covariance matrix is instead given by ((gAgA)) = + Sa) where Sa is 
some residual signal. We want to construct a test to confirm or refute the null hypothesis 
of deterministic linear bias. We aim to distinguish between the following two cases: 

• Null hypothesis Hq: ((gAgA)) = for some value of /. 

• Alternate hypothesis Hi: ((gAgA)) = Na + Sa even when / is chosen to minimize 
the contribution of Sa- 

One valid test of Hq would be a test - that is, compute a statistic assuming Hq is 
true (so = gA'^A^gA^' with respect to the parameter /, and then check to 

see whether Xmin' the value of at the minimum point, indicates a good model fit. For a 
system with N degrees of freedom where ^ 1, the distribution of can be approximated 
as a Gaussian distribution with mean A^ and variance 2A^ - here A^ = n — 1 (the number of 

cells n minus 1 fitted parameter /). Thus fxmin " ^) l^f^ 

is a measure of the number 



< 2 , i.e. 



of 'sigmas' the observed Xmin value lies away from the mean expected value if is true. 



The usual convention is to consider the model a good fit if 
the data is consistent with the null hypothesis at the 2a level. 

However, if we have some knowledge of the form that a deviation from the null hypothesis 
is likely to take (that is, a model for the residual signal Sa), we can construct a statistic that 
is more sensitive than x^ for ruling out t he null hypothesis Hn i n the e vent t hat the al t ernate 
hypothesis Hi is true. Here we follow iTegmark and Peebles (| 19981 ) and iTeemarkl (|l999l l 



aiiu g,ciiciciiiz;c Liic test statistic by constructing a quadratic test statistic q = g^Lg^ lui 

some matrix E. We want to choose E such that our test will rule out Hq with maximum 

2 



significance if Hi is true. Our generalized version of (^Xmin ~ /v 2A^, the significance 
level (i.e. the number of 'sigmas') at which we can rule out Hq^ is given by 

<l-{{<l\Ho)) 



Aq 



(6.21) 



where {AqY = {{q'^\Ho)) - {(q\Ho)) is the variance of q given that Hq is true. The desired 
matrix E maximizes 

= (6.22) 
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the expected value of v if the altern a te hy pothesis }i\ holds. The matrix E satisfying this 
requirement is derived in Tegmark ( 19991 ) to be E oc N^^SaN^^, which we insert into 



equation (j6.2ip to give our null-buster statistic: 

-TrfNT^SA) 

(6.23) 



glNX'SANX^g^-Tr(NxiSA) 



[2Tr(N^iSAN^iSA)]'^' 
This can be interpreted as the significance level (i.e. the number of 'sigmas') at which we 



can ru le out the null hypothesis of deterministic linear bias. As detailed in [Tegmark 



(|l999l l. this test assumes that the Poisson shot noise contribution can be approximated as 
Gaussian but makes no other assumptions about the probability distribution of g^. It is a 
valid test for any choice of Sa and reduces to a standard test if Sa = Na, but it rules 
out with maximum significance in the case where Ji\ is true, i.e., where Sa is the true 
residual signal. 

Using equations (16.130 . ()6.14p . and ()6.18p . the covariance matrix of gA can be written 

as 

«gAgA>> = <A if - 26reirrcl/ + bii) Sm + Na, (6.24) 

where Sm is given by equation (|6.19p . We use Sa = Sm in equation (j6.23p (note that 
z/ is independent of the normalization of S, which scales out) since deviations from linear 
deterministic bias are likely to be correlated with large-scale structure. 

To apply the null-buster test, we compute u as a function of / and then minimize it. 
If the minimum value fmin > 2, we rule out linear deterministic bias at > 2a. If the null 
hypothesis cannot be ruled out and we choose to accept it as an accurate description of the 
data, we can take rrei to be equal to unity and then use the value of / that gives as a 
measure of 6rei- 

We calculate the uncertainty on ftj-ei using two different methods. The first method 
makes use of the fact that v is generalized statistic: the uncertainty on b^ei, is given by 
the range in / that gives \/2N (z^ — t'min) ^ Ij where is the number of degrees of freedom 
(equal to the number of cells minus 1 fitted parameter). This is a generalization of the 
standard Ax^ = 1 uncertainty since v is a generalization of (^Xm m ~ ^) / V^N. The second 
method uses jackknife resampling, which is described in §6.B.1I along with a comparison of 
the two methods. We present all of our results derived from the null-buster test using the 
uncertainties from the jackknife method. 



6.3.5 Maximum Likelihood Method 



In addition to the null-buster test, we use a maximum likelihood analysis to determine the 
parameters 6rei and rj-ei. Our method is a generalization of the maximum likelihood method 
used in previous work, accounting for correlations between different cells but making a 

some what different s e t of assump t ions. 

In lBlantonI (|2000l l . Iwild et al.l (|2005l l . IConwav et al.l ^200^ ). the probability of observing 
Ni galaxies of type 1 and N2 galaxies of type 2 in a given cell is expressed as 



/oo poo 
^ J ^ Poiss [Ni,Ni{l + 5i)] 

X Poiss [N2,N2 (1 + 62)] f {61, 62, a) dJidJa, 



(6.25) 
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where f (61,62,0) is the joint probabihty distribution of 61 and 62 in one cell, a represents 
a set of parameters which depend on the biasing model, and Poiss (A^, A) = \^ /N\ is 
the Poisson probability to observe N objects given a mean value A. The likelihood function 
for n cells is given a set of observations N\ and TVg then given by 



(6.26) 



which is minimized with respect to the parameters a. This treatment makes two assump- 
tions: it neglects correlations between different cells and it assumes that the galaxy dis- 
creteness is Poissonian. These assumptions greatly simplify the computation of £, but 
are understood to be approximations to the true process. Cosmological correlations are 
known to exist on l arge scales, although their impact on count s- in-cells analyses has been 
argued to be small (Broadhurst et al. 1995 : Conwav et al. 20051 ). Semi-ana lytical modeling 



dSheth a nd Diaferioll200ll:lBerlind and Weinbergll2002l ). A'^-body simulation (ICasas-Mi randa et al 
2002 : Kravtsov et al.ll2004 ) , and smoothed particle hydrodynamic simulation ((Berlind et al 
2OO3I ) investigations suggest that the probability distribution for galaxies/halos is sub- 



Poisson ian in some regiines, and in fact nq n- Poissonian behavior is implied by observations 
as well (|Yang et al.l 1200^1 : IWild et al.ll200,5l ^. 



Dropping these two assumptions, we can write a more general expression for the likeli- 
hood function for n cells: 



.1=1 



L2 ,A^2 



(i) 



1 + 6i^^ 



2 ' 



) "2 ' 



x/(5«,...,<^i"),<5: 
xd5i'^..d5i"W(')...d4"\ 



l + 6\ 



a 



3.27) 



where Poiss (A^, A) has been replaced with a generic probability Pg {N, A, (3) for the galaxy 



distribution parameterized by some parameters /? and f (6^\ . . . ,6^\6^\ . 



,62 ,« 



is a joint probability distribution relating 61 and 62 in all cells. I n practice, this 



woul d be prohibitively difficult to calculate as it involves 2n integrations (jPodelson et al 
19971 ). and would require a reasonable parameterized form for Pg{N,X,P) as well as 



,5; 



(n) 



a 



In this chapter, we take a simpler approach and approximate the probability distri- 
bution for our data vector g to be Gaussian with the covariance matrix C as defined by 
equations (j6.12p . (j6.13p and (j6.18p . and use this to define our likelihood function in terms 
of the parameters af, b^ei, and r^ei- 
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) 




g 



(6.28) 



Note that this includes the shot noise since C = S + N, and is not precisely equivalent to 
assuming that Pg and / in equation (I6.27P are Gaussian. 

For Ti-ei values of Irreil > 1, the matrix C is singular, and thus the likelihood function 
cannot be computed. Hence this analysis method automatically incorporates the constraint 
that |rrei| < 1, which is physically expected for a cross-correlation coefficient. 

To determine the best fit values of our parameters for each pairwise comparison, we 
maximize 2 ln>C (ci^, 6rei, Trei) with respect to erf, bj-^i, and rrei- Since our method of compar- 
ing pairs of galaxy samples primarily probes the relative biasing between the two types of 
galaxies, it is not particularly sensitive to af, which represents the bias of type 1 galaxies 
relative to the dark matter power spectrum used in equation (16.19p . Thus we marginalize 
over af and calculate the uncertainty on b^ei and r^^i using the A (2 In £) = 1 contour in the 
brei-^rei plane. This procedure is discussed in more detail in §6.B.2[ 



To test the deterministic linear bias model, we apply the null-buster test described in ^6.3.41 
to the pairs of galaxy samples described in §6.3.11 For studying the luminosity-dependent 
bias, we use the galaxies in the more luminous bin as the type 1 galaxies and the dimmer 
bin as the type 2 galaxies for each pair of neighboring luminosity bins, and repeat this in 
each volume V1-V6. We do this for all galaxies and also for red and blue galaxies separately. 
For the color dependence, we use red galaxies as type 1 and blue galaxies as type 2, and 
again repeat this in each volume. To determine the scale dependence, we repeat all of these 
tests for four different values of the cell size L as described in ^6.3.21 

Is the bias linear and deterministic? 

The results are plotted in Fig. 16-61 which shows the minimum value of the null-buster test 
statistic z/jnin vs. cell size L. According to this test, deterministic linear biasing is in fact 
an excellent fit for the luminosity-dependent bias: nearly all fmin fall within |t'min| < 2, 
indicating consistency with the null hypothesis at the 2a level. (There are a few exceptions 
in the case of the red galaxies, the largest being fmin ~ 5 for the smallest cell size in V3.) 
For color-dependent bias, however, deterministic linear biasing is ruled out quite strongly, 
especially at smaller scales. 

The cases where the null hypothesis survives are quite noteworthy, since this implies 
that essentially all of the large clustering signal that is present in the data (and is visually 
apparent in Fig. I6-3P is common to the two galaxy samples and can be subtracted out. For 
example, for the V5 luminosity split at the highest angular resolution (L = 10/i^^Mpc), 
clustering signal is detected at 953fT in the faint sample (i^(/) « 953 for / = 0) and at 
255(7 in the bright sample (i£(/) ~ 255 for / = oo), yet the weighted difference of the 



6.4 Results 



6.4.1 Null-buster Results 
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two maps is consistent with mere shot noise (i^(0.88) ~ —0.63). This also shows that no 
luminosity-related systematic errors afflict the sample selection even at that low level. 

Is the bias independent of scale? 

For the luminosity-dependent bias, we use the value of / that gives fmin as a measure of 
ftrei) the relative bias between two neighboring luminosity bins. Since deterministic linear 
bias is ruled out in the case of the color-dependent bias, we instead use the value of b^ei 
from the likelihood analysis here. We find that although the value of 6rei depends on 
luminosity, it does not appear to depend strongly on scale, as can be seen in Fig. 16-71 in 
all plots the curves appear roughly horizontal. To test this 'chi-by-eye' inference of scale 
independence quantitatively, we applied a simple fit on the four data points (or three in 
the color-dependent case) in each volume using a one-parameter model: a horizontal line 
with a constant value of 6rei- For this fit we use covariance matrices derived from jackknife 
resampling, as discussed in §6.B.1[ We define this model to be a good fit if the goodness- 
of-fit value (the pro bability that a as poor as as the value calculated should occur by 
chance, as defined in lPress et ahlEooi ) exceeds 0.01. 

We find that this model of no scale dependence is a good fit for all data sets plotted 
in Fig. 16-71 We therefore find no evidence that the luminosity- or color-dependent bias 
is scale-dependent on the scales we probe here (2 — 160/i~^Mpc). This implies that re- 
cent cosmological parame t er analyses which use only measurements on scales > 60 /i^^Mpc 
(e.g. , ISanchez et aD I2OO6I : iTeemark etaP I2OO6I : ISpergel et all l2007l ^ are probably justified 



in assuming scale independence of l uminosity-dependent bias. 

In comparison to previous work dZehavi et al.ll2005l : |Li et al.ll200fil ^. it is perhaps surpris- 



ing to see as little scale dependence as we do -iLi et all (120061 ) find the luminosity-dependent 
bias to vary with scale (see their fig. 4), in contrast t o what we find here. The measurement 
of luminosity-dependent bias in IZehavi et all (|2005l ) agrees more closely with our observa- 
tion of scale independence, but their their fig. 10 indicates that we might expect to see scale 
dependence of the luminosity-dependent bias in the most luminous samples. However, we 
measure the bias in our most luminous samples (in V6) at 16 — 134/i~^Mpc, well above 



the rang e probed in IZehavi et al 



^ - jiooi), so there is no direct conflict here. Additionally, 

fig. 13 of lZehavi et al.l (|2005l ) andfig. 10 of lLi et al.1 (|2006l ) show that correlation functions of 
red and blue galaxies have significantly different slopes, implying that the color-dependent 
bias should be strongly scale-dependent on 0.1 — 10/i~^Mpc scales. However, the points 
> 1 h~^Mpc in these plots (the range comparable to the scales we probe here) do not appear 
strongly scale dependent, so our results are not inconsistent with these c orrelation function 
measurements. This interpretation is further supported by recent work ( Wang et al. 200?! ) 



that finds correlation functions for different luminosities and colors to be roughly parallel 
above ~ l/i~^Mpc. 



How bias depends on luminosity 

Our next step is to calculate the relative bias parameter 6/6* (the bias relative to 
galaxies) as a function of luminosity by combining the measured values of 6rei between the 
different pairs o f luminosity bins. Thi s function has been measured previously using SDSS 
power spectra (ITegmark et al.l 12004^ at length scal es of ^ 60 h~^ M\)c as well as SPSS 



( Zehavi et al.M200^ 



Li et al 



2006 



Wang et al.l \200i ) and 2dFGRS dNorberg et al.l l200lh 
correlation functions at length scales of ~ 1 /i~^Mpc - here we measure it at length scales 
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cell size (h ^ Mpc) cell size (h ^ Mpc) 



Figure 6-6 Null-buster results for pairwise comparisons. measures the number of sigmas 
at which deterministic linear biasing can be ruled out as a model of relative bias between 
the two samples being compared. Shaded areas indicate Iz^minl < 2, where data is consistent 
with the null hypothesis at the 2a level. Four different types of pairwise comparison are 
illustrated: (a) luminous vs. dim, (b) red vs. blue, (c) luminous red vs. dim red, and 
(d) luminous blue vs. dim blue. The different symbols denote the different comparison 
volumes V1-V6. The luminosity-dependent bias (a, c, d) is consistent with deterministic 
linear biasing but color-dependent bias (b) is not. 
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Figure 6-7 Relative bias 6rei between pairwise samples, (a) luminous vs. dim, (b) red 
vs. blue, (c) luminous r(;d vs. dim red, and (d) luminous blue vs. dim l)luc. revealing 
no significant scale dependence of luminosity- or color-dependent bias. The values 
shown for luminosity dependent splittings (a), (c), and (d) were computed with the null- 
buster analysis, those shown for the color-dependent splitting (b) were computed with the 
likelihood analysis. The different symbols denote the different comparison volumes V1-V6. 
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of - 20 /i-^Mpc. 

The bias of each luminosity bin relative to the central bin L4 is given by 



^1 A A A ^2 , , ^3 , 

— — O12O23O34, 7— — "23034, 7— — 034, 

04 04 04 

^ = \ , h = ^^ ^ = ±, (6.29) 

64 ^45^56^67 ^4 &45&56 ^4 ^45 

where denotes the measured value of 6rei between luminosity bins La and L/3 using all 
galaxies and ba denotes the bias of galaxies in luminosity bin La relative to the dark matter. 
For each pairwise comparison, we choose the value of 6rei calculated at the resolution where 
the cell size is closest to 20/i~^Mpc, as illustrated in Fig. 16-51 (Since we see no evidence 
for scale dependence of fti-ei for the luminosity-dependent bias, this choice does not strongly 
influence the results.) 

To compute the error bars on 60/64, we rewrite equation ()6.29p as a linear matrix 
equation using the logs of the bias values: 
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-1 
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/ log 61/64 \ 
log 62/64 
log 63/64 
log 65/64 
log 65/64 
V log 67/64 ) 



( log 612 \ 
log 623 
log 634 
log 645 
log 656 
V log 667 j 



(6.30) 



or Abiog 

6a/3, biog 



b 



log,rel; 



where bi, 



iog,rei is a vector of the log of our relative bias measurements 
is a vector of the log of the bias values 60/64, and A is the matrix relating them. 
We determine the covariance matrix Si-d of bj-ei (a vector of the relative bias measurements 
6q,/3) from the jackknife resampling described in §6.B.ll and then compute the covariance 
matrix Swrei of bw .el by 



^log,rel 

where Brei = diag(bi.ei). We invert equation ()6.30p to give biog: 



(6.31) 



b 



log 



log,reb 



(6.32) 



with the covariance matrix for biog given by 



■"log 



A^ yi~ A 

^ ^log.rel** 



(6.33) 



We then fit our data with the model used in Norberg et al. ( 200ll ): 6(M) /6* = ai -|- 
02 (i/-^^*), parameterized by a = (oi, 02). Here M is the central absolute magnitude of the 
bin, L is the corresponding luminosity, and = —20.83. We use a weighted least-squares 



130 




Figure 6-8 Luminosity dependence of bias for all (open circles), red (solid triangles), and 
blue (solid squares) galaxies at a cell size of ~ 20/i~^Mpc from null-buster results. The 
solid line is a mode l fit to the all-galaxy data points, the dotted line sh ows the model 
fromlTegmark et al. I (l2004bl ) . and the grey dashed line shows the model from iNorberg et al 



(|2n0lh . The lNorberg et ^1(120011 ) model has been computed using the SDSS r-band value 
of = -20.83. 



fit that is linear in the parameters (ai, 02) - that is, we solve the matrix equation 
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(6.34) 



or b = Xa, where b is a vector of the bias values ba/b^ and X is the matrix representing 
our model. We solve for a using 



a= (X^5]-iX)"'x^I]-ib. 



Here 5] is the covariance matrix of b, given by 

S = BSiogB"^, 



(6.35) 



(6.36) 



where Siog is given by equation (I6.33P and B = diag (b). This procedure gives us the best- 
fitting values for the parameters ai and 02, accounting for the correlations between the data 
points that are induced we compute the bias values b from our relative bias measurements. 
We then normalize the model such that 6 {M^) /6* = 1. 

Figure [6^ shows a plot of 6/6* vs. M: results for all galaxies are plotted with black 
open circles, our best-fitting model is shown by the solid line, the best-fitting model from 



131 




-16 -18 -20 -22 -16 -18 -20 -22 -16 -18 -20 -22 
Absolute magnitude Absolute magnitude Absolute magnitude 



Figure 6-9 Comparison to previous res u lts fo r the luminosity depe n dence of bias f or all, 
red, and blue gala xies. iNorberg et aD (|2002l l: IZehavi et all (120051'): iLi et alJ (120061 ') . and 
Wang et al. ( 200?! ) use correlation function measurements, Tegmark et al. (2004b) use the 



power spectrum, and ICroton et alJ (|2007l ) use counts in cells. To better illustrate the 
similarities and differences in the trends as a function of luminosity, we have normalized 
all measurements to match our results using the bin closest to = —20.83. The error 
bars shown are all relative: they do not include uncertainties due to the normalization. 
Numbers in parentheses denote the scale in h~^Mx>c at which t he me asurements were done. 
Also shown are theoretical models from van den Bosch et al. ( 20031 ) (we show their model 
B as a representative example) and Coorav ( 20051 ) - these are also normalized to match our 
results at M*. 



Norberg et al.l (l200lh is shown by the grey dashed line, and the best-fitting model from 
Tegmark et al. ( 2004bl ) is shown by the dotted line. The error bars represent the diagonal 
elements of S from equation (j6.36p . Our model, with (a i,a^) = (0.862, 0.138), agrees 
extremely well with the model from Norberg et al. ( 200ll ). with (01,02) = (0.85, 0.15). 
This agreement is quite remarkable since we use data from a different survey and analyze 
it with a completely different technique. 



A comparison of our results with previous measurements is shown in Fig. 16-91 in the left 



pane l . In order to comp are our SDSS results with results from 2dFGRS ([Norberg et al 
2OO2I : ICroton et al.ll2007l ). we have added a constant factor o f —1.13 to their quote d values 



for Mfjj — 5 logio ^ order to line up the value of used in iNorberg et all (|2002l l (M, 



51ogio h = —19.7) with the value used here (Mo: 



-20.83). Note that this is necessarily 



a rough correction since the magnitude in the different bands varies depending on the 
spectrum of each galaxy, but this method provides a reasonable means of comparing the 
different results. This plot shows excellent agreement over a wide range of scales, lending 
further support to our conclusion that the luminosity-dependent bias is independent of 
scale. 



We also use equation (j6.32p to calculate 6/6* vs. M for red and blue galaxies separately. 
To plot the points for red, blue, and all galaxies on the same 6/6* vs. M plot, we need 
to determine their relative normalizations. Applying equation (j6.32p to the red and blue 
galaxies gives 6Q^red/64,red and 6a,biue/64,biue, so to normalize the red-galaxy and blue-galaxy 
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data points to the all-galaxy data points in Fig. 16-81 we need to calculate 

^a,red ^4, all ^4, red ^«,red 



^*,all 



&*,all ^'4,all ^4,red 



(6.37) 



and 



^a,blue _ 04,all 04,red 04,blue Og.blue 
^*,all K,all &4,all &4,red ^4,blue 



(6.38) 



The factor 64,aii/^*,aii is simply the normalization factor chosen for the above model to 
give b{M^) /b^: = 1. To determine ?)4,red/^4,aib we use best-fitting values of from the 
likelihood analysis described in ^6.3.5l at the resolution with cell sizes closest to 20/i~^Mpc: 
o"! from the comparison of dimmer and more luminous galaxies in V3 gives 64_aii) and 
similarly cxi from the comparison of blue and red galaxies in L4 gives 64, red) so 



b. 



4,red 
^4,all 



2 \ 1/2 

"^l,redvs. blueL4 \ 



l,lum vs. dim V3 ^ 



(6.39) 



The blue points are then normalized relative to the red points using 64,biue/^4,red equal to the 
measured value of ^i-ei from the likelihood comparison of blue and red galaxies in L4. Thus 
the shapes of the red and blue curves are determined using the luminosity-dependent bias 
from the null-buster analysis, but their normalization uses information from the likelihood 
analysis as well. 

Splitting the luminosity dependence of the bias by color reveals some interesting fea- 
tures. The bias of the blue galaxies shows only a weak dependence on luminosity, and 
both luminous (M ~ —22) and dim (M ~ —17) red galaxies have slightly higher bias than 
moderately bright (M ~ —20 ~ M*) red galaxies. The previously observed luminosity 
dependence of bias, with a weak dependence dimmer than L^, and a strong increase above 
L*, is thus quite sensitive to the color selection: the lower luminosity bins contain mostly 
blue galaxies and thus show weak luminosity dependence, whereas the more luminous bins 
are dominated by red galaxies which drive the observed trend of more luminous galaxies 
being more strongly biased. It is instructive to compare these results with the mean local 
overdensity in color-magnitude space, as in fig. 2 of Blanton et al. (2005a). Although our 
bias measurements are necessarily much coarser, it can be s een that the bias is stro ngest 
where the overdensity is largest, as has been seen previously ( Abbas and Sheth 20061 ). 

Comparisons of our results to other measurements of luminosity-dependent bias for red 
and blue galaxies are shown in Fig. 16-91 in the middle and right panels. Indications of the 
differing trends for red and blue galaxies have been observed in previous wo r k: an early 
hint of the upt urn in the bias fo r dim red galaxies was seen in iNorberg et al] ( 20021 ) . and 
recent results ( Wang et al. 20071 ) also indicate higher bias for dim red galaxies at scales 
> l/i~^Mpc . How ever, there is some inconsistency between these results compared to 
Zehavi et al.] (|2005l l and lLi et al.l toodi ) regarding the dim red galaxies: they find that dim 
red galaxies exhibit the strongest clustering on scales < 1 h~^Mpc and luminous red galax- 
ies exhibit the strongest clustering on larger sca l es, as can be seen fr om the g r een p oints in 
Fig. [EH This is shown in fig. 14 of lZehavi etdH (|2005l l and fig. 11 of lLi et al.l (jiooi). How- 
ever, we find the dim red galaxies to have higher bias than L* red galaxies at all the scales 
we probe (2—40 /i~^Mpc in this case). This uptu rn of the bias for dim red gal axies is present 
in the halo model-based theoretical curves from van den Bosch et al.l (|2003l ). although not 
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in the theoretical curves from Cooravl Note also that Ivan den Bosch et alJ (|2003l l 

use the data from lNorberg et alJ (2002) to constrain their models so the agreement between 
the theory and data should be interpreted with some caution. 



Recent measurements of higher order clustering statistics (jCroton et alJ 120071 ) find the 
same trends in the clustering strengths of red and blue galaxies, although they indicate that 
their linear bias measurement (which should be comparable to ours) shows the opposite 
trends - little luminosity dependence for red galaxies and a slight monotonic increase for 
blue galaxies. However, their luminosity range is much narrower than ours so the trends are 
less clear, and placing their data points on Fig. 16-91 shows that they are in good agreement 
with our results. 

Previous studies (|Norberg et al.1 12OO2I : iLi et al.l I2OO6I : IWang et all booj ) have also re- 
ported a somewhat stronger luminosity depe ndence of blue galaxy clus tering than we have 



measured here. As can be seen in Fig. 16-9^ iNorberg et al.l (120021') and IWang et al.l (120071 ) 



measure slightly higher bias for luminous blue galaxies, and (|Li et al.ll2006l ) measure slightly 
lower bias for dim blue galaxies. Although the quantitative disagreement is fairly small, the 
qualitative trends of the previous studies imply that the bias of blue galaxies increases with 
luminosity, as opposed to our measurement which indicates a lack of luminosity dependence. 



6.4.2 Likelihood Results 

To study the luminosity dependence, color dependence and stochasticity of bias in more 
detail, we also apply the maximum likelihood method described in ^6.3.5l to all of the same 
pairs of samples used in the null-buster test. Due to constraints on computing power and 
memory, we perform these calculations for only three values of the cell size L rather than 
four, dropping the highest resolution (smallest cell size) shown in Fig. 16-41 The likelihood 
analysis makes a few additional assumptions, but provides a valuable cross-check and also 
a measurement of the parameter rrei which encodes the stochasticity and non-linearity of 
the relative bias. 

For each pair of samples, the likelihood function given in equation (j6.28p is maximized 
with respect to the parameters af, 6rei) and r^ei and marginalized over erf to determine the 
best-fitting values of b^ei and rj-ei, with uncertainties defined by the A (21ni2) = 1 contour 
in the 6rei-rrei plane. As we discuss in ^6.B.2l the values of bj-^i found here are consistent 
with those determined using the null-buster test. 

Figure 16-101 shows the best-fitting values of function of cell size L. For the 

comparisons between neighboring luminosity bins, the results are consistent with r^-d = 1. 
On the other hand, the comparisons between red and blue galaxies give rj-ei < 1, with 
smaller cell sizes L giving smaller values of rrei. This confirms the null-buster result that 
the luminosity-dependent bias can be accurately modeled using simple deterministic linear 
bias but color-dependent bias demands a more complicated model. Also, r^^i for the color- 
dependent bias is seen to depend on scale but not strongly on luminosity. In contrast, ftj-gj 
(both in the null-buster and likelihood analyses) depends on luminosity but not on scale. 

To summarize, we find that the simple, deterministic model is a good fit for the luminosity- 
dependent bias, but the color-dependent bias shows evidence for stochasticity and/or non- 
linearity which increases in strength towards smaller scales. These results are consis- 
tent with previous detections of stochasticity /non-linearity in s pectral-type-depende nt bias 
(ITegmark and Bromlevl[l999l : iBlantonllioOol : IConwav et al.boOfil ) , and also agree with (jWild et al. 
2OO5I ) which measures significant stochasticity between galaxies of different color or spectral 
type, but not between galaxies of different luminosities. 
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Figure 6-10 The best-fitting values of the relative cross-correlation coefficient r^c\ between 
pairwise samples. Four different types of pairwise comparison are illustrated: (a) luminous 
vs. dim, (b) red vs. blue, (c) luminous red vs. dim red, and (d) luminous blue vs. dim 
blue. The different symbols denote the different comparison volumes V1-V6. 
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Figure 6-11 Comparison of relative cross-correlation coefficient rrei 
galaxies as measured with different tecliniques. The points from IZehavi et al 



between red and blue 
(|2005l l 



are 



extracted from cross-co rrelation measurem ents between red and blue galaxies in SDSS with 
Mo.ir < —21, and the I Wild et alj (j2005l ) points are f r om a. counts- in-cells analysis using 



all 2dFGRS galaxies. Our results and the I Wang et alJ (120071 ) results (also fror n SDSS) are 
separa ted by luminosity - symbols are the same as in Fig. 16-101 and for the I Wang et al. 
(|2007l ) results, open triangles denote their dimmest bin (—19 < Moi.j. < —18) and open 
sq uares denote their most luminous bin (—23 < Mol^ < —21.5). The length scales used 
Wang et al.l (j2007l ) are averages over small scales (0.16 — 0.98 /i ^Mpc) and large scales 



m 
(0.98 - 
clarity. 



9.8 h ^Mpc) - points here are shown in the middle of these ranges and offset for 



We compare our results for r^^i for red and blue galaxies to previous re s ults in Fig. 16-111 
This shows good agreement between our results and those of Wild et al. I (|200,5l ) (n in from 
their fig. 11), implying that these results are quite robust since our analysis uses a different 
data set, employs different methods, and makes different assumptions. 

For the results f r om cr oss-correlation measurements, however, the agreement is not as 
clear. IZehavi et al.l (120051 ) find that the cross-correlation between red and blue galaxies 
(their fig. 24), indicates that rj-ei is consistent with 1 on scales > 1 /i^^Mpc. However, it 
is not clear that this result disagrees with ours, as their result is for luminous galaxies 
{Mo.ij. < —21) and and we do not see a strong indication of r^-d < 1 for our V6 sample 
(23 < Mo.i^ < -21). 

More recent cross-correlation measurements (jWang et al.l 120071 ) do find evidence for 
stochasticity /non-linearity between red and blue galaxies at scales < 1/i^^Mpc and also 
show an indication that dimmer galaxies have slightly lower values of rjei. Note also that 
the met hod of calcul a ting r rpi b y taking ratios of cross- and auto-correlation functions as 
used for lZehavi et"aD (|2005l ) and I Wang et~all jioO^) does not automatically incorporate the 
constraint that |rrei| < 1 as our analysis does, so their error bars are allowed to extend 
above Trei = 1 in Fig. 16-111 

Overall, the counts-in-cells measurements (this chapter. Wild et al. 20051 ) show stronger 
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evidenc e for stochasticity/ non-linearity at lar ger scales than the cross-correlation measure- 
ments ( Zehavi et al. 20051 : Wang et al. 2007 ). indicating either that there might be some 
slight systematic variation between the two methods or that the counts-in-cells method is 
more sensitive to these effects. 



6.5 Conclusions 



To shed further light on how galaxies trace matter, we have quantified how different types 
of galaxies trace each other. We have analyzed the relative bias between pairs of volume- 
limited galaxy samples of different luminosities and colors using counts-in-cells at varying 
length scales. This method is most sensitive to length scales between those probed by 
correlation function and power spectrum methods, and makes point-by-point comparisons of 
the density fields rather than using ratios of moments, thereby eliminating sample variance 
and obtaining a local rather than global measure of the bias. We applied a null-buster test 
on each pair of subsamples to determine if the relative bias was consistent with deterministic 
linear biasing, and we also performed a maximum-likelihood analysis to find the best-fitting 
parameters for a simple stochastic biasing model. 

6.5.1 Biasing results 

Our primary results are: 



1. The luminosity-dependent bias for red galaxies is significantly different from that of 
blue galaxies: the bias of blue galaxies shows only a weak dependence on luminosity, 
whereas both luminous and dim red galaxies have higher bias than moderately bright 
(L^) red galaxies. 

2. Both of our analysis methods indicate that the simple, deterministic model is a 
good fit for the luminosity-dependent bias, but that the color-dependent bias is more 
complicated, showing strong evidence for stochasticity and/or non-linearity on scales 
< lOh-^ Mpc. 



3. The luminosity-dependent bias is consistent with being scale-independent over the 
range of scales probed here (2 — 160/i~^Mpc). The color-dependent bias depends on 
luminosity but not on scale, while the cross-correlation coefficient rrei depends on scale 
but not strongly on luminosity, giving smaller Tfei values at smaller scales. 



These results are encouraging from the perspective of using galaxy clustering to measure 
cosmological parameters: simple scale-independent linear biasing appears to be a good ap- 
proximation o i i the ^. 60/t~^Mpc scales u s ed in many recent cosmological studies (e.g., 
Sanchez et all (|2006l ): lTegmark et all (|2006l ^: ISDergel et al.l tooii )). However, further quan- 



tification of small residual effects will be needed to do full justice to the precision of next- 
generation data sets on the horizon. Moreover, our results regarding color sensitivity sug- 
gest that more detailed bias studies are worthwhile for luminous red galaxies, which have 
emerged a powerful cosm ological probe because of the i r visibility at large di stances and near- 
optimal number density ( Eisenstein et al. 2001 . 2005 : Tegmark et al. 20061 ). since color cuts 
are involved in their selection. 
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6.5.2 Implications for galaxy formation 



What can these results tel l us about ga l axy f ormation in the context of the halo model? 
First of all, as discussed in IZehavi et all (|2005l ^. the large bias of the faint red galaxies can 
be explained by the fact that such galaxies tend to be satellites in high mass halos, which 
are more strongly clustered than low mass halos. Previous studies have found that central 
galaxies in low-mass halos are preferentially blue, central galaxies in high mass halos tend 
to be red, and that the luminosity of the c entral galaxy is strongly correlated with the halo 
mass dYang et alJ 1200,4 Izheng et al.ll200,^ ). Our observed lack of luminosity dependence of 
the bias for blue galaxies would then be a reflection of the correlation between luminosity 
and halo mass being weaker for blue galaxies than for red ones. Additional work is needed to 
study this quantitatively and compare it with theoretical predictions from galaxy formation 
models. 

The detection of stochasticity between red and blue galaxies may imply that red and blue 
galax ies tend to live in different halos - a study of galaxy groups in SDSS (jWeinmann et al 



20061 ) recently presented e vidence supp o rting this, but this is at odds with the cross- 
correlation measurement in IZehavi et al.] ^200^ ). which implies that blue and red galaxies 
are well- mixed within halos. The fact that the stochasticity is strongest at small scales 
suggests that this effect is due to the 1-halo term, i.e., arising from pairs of galaxies in the 
same halo, although some amount of stochasticity persists even for large scales. However, 
the halo model implications for stochasticity have not been well-studied to date. 

In summary, our results on galaxy biasing and future work along these lines should be 
able to deepen our understanding of both cosmology (by quantifying systematic uncertain- 
ties) and galaxy formation. 
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Figure 6-12 Null-buster results for randomly split samples. 

6. A Consistency Checks 

6.A.1 Alternate null-buster analyses 

In order to test the robustness of our results against various systematic effects, we have 
repeated the null-buster analysis with four different modifications: splitting the galaxy 
samples randomly, offsetting the pixel positions, using galaxy positions without applying 
the finger-of-god compression algorithm, and ignoring the cosmological correlations between 
neighboring cells. 



Randomly split samples 

The null-buster test assumes that the Poissonian shot noise for each type of galaxy in each 
pixel is uncorrelated - i.e., that the matrix N in equation (|6.13|) is diagonal - and that this 
shot noise can be approximated as Gaussian. To test the impact of these assumptions on 
our results, we repeated the null-buster analysis using randomly split galaxy samples rather 
than splitting by luminosity or color. For each volume VI- V6, we created two samples 
by generating a uniformly distributed random number for each galaxy and assigning it to 
sample 1 for numbers > 0.5 and sample 2 otherwise. 

If the null-buster test is accurate, we expect the pairwise comparison for the randomly 
split samples to be consistent with deterministic linear bias with ^j-ei = 1- The results are 
shown in Fig. 16-12! ~ we find that, indeed, deterministic linear bias is not ruled out, with 
nearly all of the Umin points falling within ±2. Furthermore, the measured values of b^ei are 
seen to be consistent with 1. Thus, we detect no systematic effects due to the null-buster 
assumptions. 

Offset pixel positions 

To test if our results are stable against the pixelization chosen, particularly at large scales 
where we have a small number of cells, we shifted the locations on the sky of the angular 
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pixels defining the cells by half a pixel width in declination. Applying the null-buster 
analysis to the offset pixels reveals no significant differences from the original analysis: the 
luminosity-dependent comparisons for all, red, and blue galaxies are still consistent with 
deterministic linear bias, and the color-dependent comparison still shows strong evidence 
for stochasticity and/or nonlinearity, especially at smaller scales. 

For the all-, red-, and blue-galaxy, luminosity-dependent comparisons, we also compared 
the measured values of ^i-ei from the offset analysis with the original analysis: we took the 
difference between the two measured values in each volume at each resolution and divided 
this by the larger of the error bars on the two analyses to determine the number of sigmas 
by which the two analyses differ. In order to be conservative, we did not add the error 
bars from the two analyses in quadrature, since this would overestimate the error on the 
difference if they are correlated and this would make our test less robust. Note that a 
fully proper treatment would necessitate accounting for the correlations between the errors 
from each analysis, which we have not done - the discussions in the this and the following 
sections are meant only to serve as a crude reality check. 

For the error bars on the original analysis, we used the jackknife uncertainties described 
in ^6.B.ll and for the offset analysis error bars we use the generalized uncertainties 
described in §6.3.41 computed from the offset results. (This is because we did not perform 
jackknife resampling for the offset case or the other modified analyses.) 

The results show good agreement: out of a total of 72 measured b^ei values, only 4 differ 
by more than 2a (all galaxies in V3 and V5 at the second-smallest cell size, at 2.6a and 3.2a 
respectively, and the red galaxies in V3 and V5 at the second-smallest cell size, at 2.2cTand 
3.3a respectively). As a rough test for systematic trends, we also counted the number of 
measurements for which the measured value of bj-e\ is larger in each analysis. We found that 
in 38 cases the value from the offset analysis was larger, and in 34 cases the value from the 
original analysis was larger, indicating no systematic trends in the deviations. 



No finger-of-god compression 



Our analysis used the finger-of-god compression algorithm from Tegmark et al. ( 2004bl ) 
with a threshold density of 5c = 200. This gives a first-order correction for redshift space 
distortions , but it cornplicat es comparisons to other analyses which wo rk purely in redshif t 
space (e.g. Iwild et ahlEioi ^ or use projected correlation functions (e.g. IZehavi et al.ll2005l ). 
particularly at small scales (< 10/i~^Mpc) where the effects of virialized galaxy clusters 
could be significant. To test the sensitivity of our results to this correction, we repeated 
the null-buster analysis with no finger-of-god compression. 

The results show excellent agreement with the original analysis - the smallest-scale 
measurements for the color-dependent comparison only rule out deterministic linear bias at 
30 sigma rather than 40, but the conclusions remain the same. Additionally, we compared 
the measured ^j-ei values to the original analysis as in ^6. A. II and find all 72 measurements 
to be within 2a. In 25 cases, the analysis without finger-of-god compression gave a larger 
bj-ei value, and in 47 cases the original analysis gave a larger value. This indicates that there 
might be a very slight tendency to underestimate 6rei if fingers-of-god are not accounted 
for, but the effect is quite small and well within our error bars. Thus, the finger-of-god 
compression has no substantial impact on our results. 
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Uncorrelated signal matrix 



The null-buster test requires a choice of residual signal matrix S a ~ our analysis uses a signal 
matrix derived from the matter power spectrum, thus accounting for cosmological correla- 
tions between neighboring cells. However, t hese correlations are comrnonly assumed to be 
negligible in other counts- in-cells analyses ( Blanton 2000l : Wild et al. 2005 : Conway et al 



I2OO5). To test the sensitivity to the choice of Sa, we repeated the analysis using Sa equal 
to the identity matrix. 

Again, we find the results to agree well with the original analysis and lead to the same 
conclusions. When comparing the measured 6rei values to the original analysis as in §6.A.H 
we find only 4 out of 72 points differing by more than 2a (all galaxies in V5 at the smallest 
cell size, at —2. 60", red galaxies in V4 at the second-smallest and smallest cell size, at 2.2a 
and 2.8a respectively, and blue galaxies in V4 at the smallest cell size, at — 3.1cj). In 39 
cases the value of b^ei is larger with the uncorrelated signal matrix, and in 33 cases b^ei is 
larger in the original analysis, indicating no strong systematic effects. Thus we expect our 
results to be directly comparable to other counts-in-cells analyses done without accounting 
for cosmological correlations. 



6.B Uncertainty Calculations 

6.B.1 Jackknife uncertainties for null-buster analysis 

We use jackknife resampling to calculate the uncertainties for the null-buster analysis. The 
concept is as follows: divide area covered on the sky into N spatially contiguous regions, and 
then repeat the analysis times, omitting each of the A^ regions in turn. The covariance 
matrix for the measured parameters is then estimated by 

^ = ^ E - ^Ll) {Kel,k - Ke^) (6-40) 

k=l 

where superscripts i and j denote measurements of b^ei in different volumes and at different 
scales, 6*gj ^ denotes the value of 6*^; with the kth jackknife region omitted, and is the 
average over all A^ values of 6*^^ ^ . 

For our analysis, we use the 15 pixels at our lowest resolution (upper left panel in 
Fig. l6-4p as the jackknife regions. However, since we use a looser completeness cut at the 
lowest resolution, two of these pixels cover an area that is not used at higher resolutions. 
Thus we chose not to include these two pixels in our jackknifes since they do not omit much 
(or any) area at the higher resolutions. Thus our jackknife resampling has A^ = 13. This 
technique allows us to estimate the uncertainties on all of our 6rei measurements as well as 
the covariance matrix quantifying the correlations between them. We use these covariances 
in the model-fitting done in §6.4.11 and §6.4.11 

Figure [6- 1 3 1 shows the uncertainties on ^rei as calculated from jackknife resampling com- 
pared to those calculated with the generalized method described in ^6.3.41 Overall, the 
two methods agree well, but the jackknife method gives larger uncertainties at the smallest 
scales and in volume VI. The reason for the large jackknife uncertainties in volume VI 
is because it is significantly smaller than the other volumes, and it is small enough that 
omitting a cell containing just one large cluster can have a substantial effect on the mea- 
sured value of 6rei- Thus the large uncertainties in VI reflect the effects of sampling a small 
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Figure 6-13 Comparison of two methods for calculating uncertainties on fti-ei from the null- 
buster analysis: jackknife resampling (red) and the generalized method (black). Also 
shown are the results for 6rei from the likelihood analysis (grey). 
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Figure 6-14 Typical contour plots of A (21n£) for volume V5 for three different resolutions 
corresponding (from left- to right-hand-side) to cell sizes of 89, 39, and 21 /i^^Mpc. Blue 
contours denote the 1, 2 and 3a two-dimensional confidence regions, and black contours de- 
note the lo" one-dimensional confidence region used for computing error bars on ftj-ei and ri-ei- 
The red contours denote the error ellipse calculated from the second-order approximation 
to 21n£ at the best-fit point, marked with a x. 



volume. Since there are so few dim red galaxies, these effects are particularly egregious for 
the measurements of luminosity-dependent bias of red galaxies in VI. Thus, based on the 
jackknife results, we elected to not use VI in our analysis of the red galaxies. 



6.B.2 Likelihood uncertainties 
Likelihood contours 

As described in §6.3.51 we calculate the uncertainty on 6rei and rj-ei for the likelihood method 
using the A (2 In £) = 1 contour in the ftrel-'^rel plane after marginalizing over cr^. This means 
that for each comparison volume and at each resolution, we calculate C from equation (j6.28p 
over a grid of 6rei and rrei values and maximize 2lnC with respect to at each grid point. 
This gives us a 2-dimensional likelihood function, which we then maximize to find the 
best-fit values for ^^ei and rj-ei- Uncertainties are calculated using the function 

A (2 In £ (6rci rrci)) ^ 2 In £ {b^fy^D - 2 ln£ (6rei r^ei) • (6.41) 

Typical contour plots of this function for volume V5 at each of the three cell sizes used are 
shown in Fig. I6-14[ 

We d efine 1- and 2-dini ensional confidence regions using the standard procedures de- 
tailed in Press et al. (|l992l ). using A(21n£) as an equivalent to Ax^: the la (68.3%) 



1-dimensional confidence region is given by A(21n£) = 1, so we define our error bars on 
6rel and ri-ei by projecting the A (21n£) = 1 contour (shown in black in Fig. I6-14P onto the 
6rei and Tj-ei axes. For illustrative purposes we also show the la, 2a, and 3a (68.3%, 95.4%, 
and 99.73%) 2-dimensional confidence regions in these plots, given by A(21n£) = 2.30, 
6.17, and 11.8 respectively. 

To check the goodness of fit, we also compute an effective value of x^'- 

Xeff = -21n£-ln|C| -2nln(27r), (6.42) 
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where C is given by equation (j6.12p and n is the number of cells. If our model is a good fit, 
the value of xis hest fit parameter values should be close to the number of degrees of 

freedom, given by dof = 2n — 2 (2n data points for type 1 and 2 galaxies in each cell minus 2 
parameters 6rei and rrei). We calculated xls/doi for each volume and resolution and found 
they all lie quite close to 1, ranging from a minimum value of 0.678 to a maximum value of 
1.11. Thus this test indicates our model is a good fit. 

The uncertainties on ^j-ei and r^^i could perhaps be calculated more accurately using 
jackknife resampling as we did for the null-buster case; however, repeating the analysis for 
each jackknife sample is computationally prohibitive since performing all the calculations 
for just one likelihood analysis took several months of CPU time. 



brei-Vrei covariance matrices 



Alternatively, we can calculate the uncertainties using the parameter covariance matrix at 
the best fit parameter values, as is commonly done in analyses. The Hessian matrix of 
second derivatives is given by 

( d^(21n£) d^(21n£) \ 
^ 2(2f:S (6-43) 

and the parameter covariance matrix is given by 

Cparam = 2H-^ = I ) . (6.44) 

Thus the uncertainties are given by ^ and a^^^^ with this method. This is equivalent to 
approximating the likelihood function 21n£ with its second-order Taylor series about the 
best-fit point, and it defines an error ellipse that approximates the A (21n£) = 1 contour. 
These error ellipses are shown in Fig. 16-141 in red, and are seen to be in close agreement 
with the true A (21n£) = 1 contours. 

This method also allows us to measure the correlation between 5rei and r^ei by calculating 
the correlation coefficient, given by 

2 

R = ^'''^"''"^"^ ^ (6.45) 

^2 0-2 ^ ' 

fcrcl ''rol 

R will fall between -1 (perfectly anti-correlated) and 1 (perfectly correlated). Effectively 
this measures the tilt of the error ellipse in the 6rei-^rei plane. Overall we find the values 
of R to be quite small - typically \R\ ~ 0.05 - indicating no large correlations between 
6rei and rrd- Out of the 72 points we calculate, only 6 have \R\ > 0.3. The cases with 
the largest R values are for blue galaxies in volume V6, where the uncertainties are quite 
large due to the small number of bright blue galaxies and the error ellipses are not good 
approximations to the likelihood contours anyway - thus these few cases with large R are 
not overly concerning. 
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Comparison with null-buster results 

Finally, we compare the results for bj-d from the likelihood method to the results from the 
null-buster analysis in Fig. 16-131 with the likelihood points shown in grey. As can be seen 
in this plot, the likelihood and null-buster values for ^j-ei agree within the uncertainties, 
even for the color-dependent bias where the null-buster values are not necessarily accurate 
since deterministic linear bias is ruled out. Thus our two analysis methods are in excellent 
agreement with each other. 
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Chapter 7 

Conclusions 



We set out in this thesis to explore the connections between particle physics and astrophysics 
by using the tools of particle physics to do astronomy, a.k.a. "Astrophysics Underground," 
and by using the tools of astronomy to do particle physics, a.k.a. "Particle Physics in the 
Sky." What have we accomplished in these endeavors? 

In "Astrophysics Underground," we performed a search for extremely high-energy upward- 
going muons in the Super-K neutrino detector and set an upper limit on the flux of neutrinos 
coming from astrophysical sources such as AGNs. Have we really learned anything about 
astrophysics here? Yes, albeit rather indirectly - our upper limit on the astrophysical 
neutrino flux tells us that AGNs are not producing an unexpectedly large amount of neu- 
trinos. However, much more time and energy was devoted to thinking about the particle 
physics technical issues than about the primary astrophysical science goal. We developed a 
way to use previously ignored information from the Super-K detector, studied how muons 
propagate through rock and water, and did a detailed analysis of various particle physics 
phenomena contributing to the background flux. Ultimately, most of what we learned was 
really about particle physics. 

This same effect was even more true in "Particle Physics in the Sky." Here we developed 
a powerful software tool to manage angular masks of the next generation of galaxy surveys 
and analyzed the details of relative bias between galaxies of different luminosities and colors. 
These projects contribute to cosmology by providing techniques to better use galaxy surveys 
as a cosmological tool. In particular, gaining a better understanding of bias is key for 
reducing the systematic uncertainties inherent in galaxy survey cosmology. Have we learned 
any particle physics? Again, yes, though again rather indirectly, in that we have contributed 
to the overall cosmological programs to understand dark matter and dark energy. What 
we really learned about is the physics of galaxies - our results on how different types 
of galaxies cluster differently most directly shed light on the complex physical processes 
governing galaxy formation and evolution. It may appear overly optimistic to claim that 
we have really done any particle physics here, but it is important to remember that our 
understanding of cosmology from galaxy surveys is only as good as our understanding of 
bias, and if we hope to gain further insight about particle physics from the next generation 
of cosmological measurements, it is essential that we understand the physics of galaxies at 
a deeper level. 

One particularly interesting direction for future research is to use relative bias measure- 
ments as a guide for refining halo model methods for putting galaxies in dark matter halos, 
as discussed in §4.1.51 This can help us gain a better understanding of the physics that 
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links galaxy properties to halo mass, and provide a better way of modeling of complicated 
bias. Furthermore, we can combined this with advances in nonlinear modeling of structure 
formation in order to further refine galaxy surveys as a tool for doing cosmology. This makes 
the biggest difference at length scales slightly smaller than where the power spectrum is typ- 
ically measured, a regime which is particularly important for setting cosmological bounds 
on the neutrino mass. Such investigations will also help reduce systematic uncertainties on 
measuring the baryon oscillation length scale as a standard ruler for studying dark energy. 
Thus, even if we have not addressed them directly, the particle physics science goals are 
clearly in sight. 

It is perhaps a natural outcome of studying the science questions of one field with the 
techniques of another that one gets more immersed in the area of the technique. This is 
borne out in the traditional physics department divisions as well - the work on neutrino 
astrophysics presented here was done within a particle physics division, and in general is 
done by scientists who call themselves particle physicists rather than astronomers. Likewise, 
the work on cosmology was done within an astrophysics division and is more often grouped 
with astronomy, not high energy physics. In many ways this makes sense - after all, if a 
scientist at such an interface spends most of their time thinking about the physics of the 
techniques they are using, it is perfectly logical to surround oneself with experts in that 
area. However, more direct interaction between scientists on either side of the fence can 
potentially be extremely beneficial, as is happening at new interdisciplinary research centers 
across the world. 

In conclusion, we have hopefully provided valuable tools for upcoming experiments in 
both neutrino astronomy and cosmology. These lines of research represent areas where the 
traditionally separate fields of particle physics and astrophysics overlap and have produced 
much fruitful interaction between physicists exploring the largest and smallest scales in our 
universe. The work in this thesis makes contributions to the exciting research in progress 
at the interface of these two rich fields. 
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